diff --git a/INIT/ftINIT.m b/INIT/ftINIT.m index e9d6a9b4..29654b29 100644 --- a/INIT/ftINIT.m +++ b/INIT/ftINIT.m @@ -54,9 +54,7 @@ % useScoresForTasks true if the calculated reaction scored should be % used as weights when fitting to tasks (optional, default % true) -% paramsFT parameter structure as used by getMILPParams. This -% is for the fitTasks step. For the INIT algorithm, -% see params (optional, default []) +% paramsFT *obsolete option* % verbose if true, the MILP progression will be shown. % (optional, default false) % diff --git a/INIT/ftINITFillGaps.m b/INIT/ftINITFillGaps.m index 5673c3e6..6b1aa471 100644 --- a/INIT/ftINITFillGaps.m +++ b/INIT/ftINITFillGaps.m @@ -23,7 +23,7 @@ % reference tModel. % The solver will try to maximize the sum of the % scores for the included reactions -% params parameter structure as used by getMILPParams +% params *obsolete option* % verbose if true, the MILP progression will be shown. % % addedRxns the rxns added diff --git a/INIT/ftINITFillGapsForAllTasks.m b/INIT/ftINITFillGapsForAllTasks.m index 0fdbe898..fbee4782 100644 --- a/INIT/ftINITFillGapsForAllTasks.m +++ b/INIT/ftINITFillGapsForAllTasks.m @@ -18,7 +18,7 @@ % reactions (optional, default is -1 for all reactions) % taskStructure structure with the tasks, as from parseTaskList. If % this is supplied then inputFile is ignored -% params parameter structure as used by getMILPParams +% params *obsolete option* % verbose if true, the MILP progression will be shown. % % diff --git a/INIT/ftINITFillGapsMILP.m b/INIT/ftINITFillGapsMILP.m index 45fb5a9f..cc66f33c 100644 --- a/INIT/ftINITFillGapsMILP.m +++ b/INIT/ftINITFillGapsMILP.m @@ -12,7 +12,7 @@ % with the same number of elements as reactions in the model, % of a vector of indexes for the reactions that should be % minimized (optional, default model.rxns) -% params parameter structure as used by getMILPParams (optional) +% params *obsolete option* % scores vector of weights for the reactions. Negative scores % should not have flux. Positive scores are not possible in this % implementation, and they are changed to max(scores(scores<0)). diff --git a/INIT/getINITModel.m b/INIT/getINITModel.m index e48bf1a2..37c16e9f 100644 --- a/INIT/getINITModel.m +++ b/INIT/getINITModel.m @@ -1,10 +1,10 @@ -function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure, params, paramsFT) +function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure) % getINITModel_legacy % Generates a model using the INIT algorithm, based on proteomics and/or % transcriptomics and/or metabolomics and/or metabolic tasks. This is the original % implementation of tINIT, which is replaced by ftINIT. % -% Input: +% Input: % refModel a model structure. The model should be in the % closed form (no exchange reactions open). Import % using import(filename,false). If the model is not @@ -53,15 +53,8 @@ % as an alternative way to define tasks when Excel % sheets are not suitable. Overrides taskFile (optional, % default []) -% params parameter structure as used by getMILPParams. This is -% for the INIT algorithm. For the the MILP problems -% solved to fit tasks, see paramsFT (optional, default []) -% paramsFT parameter structure as used by getMILPParams. This is -% for the fitTasks step. For the INIT algorithm, see -% params (optional, default []) % -% -% Output: +% Output: % model the resulting model structure % metProduction array that indicates which of the % metabolites in metabolomicsData that could be diff --git a/INIT/runINIT.m b/INIT/runINIT.m index 33ddaac0..ca064f50 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,params) +function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops) % 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,8 +38,6 @@ % problem significantly more computationally intensive to % solve (two more integer constraints per reversible reaction) % (optional, default false) -% params parameter structure as used by getMILPParams (optional, -% default []) % % outModel the resulting model structure % deletedRxns reactions which were deleted by the algorithm diff --git a/core/consumeSomething.m b/core/consumeSomething.m index a9c7f1ee..5ec896cb 100755 --- a/core/consumeSomething.m +++ b/core/consumeSomething.m @@ -19,7 +19,7 @@ % fluxes instead of the sum. Slower, but can be % used if the sum gives too many fluxes (optional, default % false) -% params parameter structure as used by getMILPParams (optional) +% params *obsolete option* % ignoreIntBounds true if internal bounds (including reversibility) % should be ignored. Exchange reactions are not affected. % This can be used to find unbalanced solutions which are diff --git a/core/fillGaps.m b/core/fillGaps.m index bf19411b..da0d81c9 100755 --- a/core/fillGaps.m +++ b/core/fillGaps.m @@ -1,4 +1,4 @@ -function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params) +function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores) % fillGaps % Uses template model(s) to fill gaps in a model % @@ -24,7 +24,6 @@ % The solver will try to maximize the sum of the % scores for the included reactions (optional, default % is -1 for all reactions) -% params parameter structure as used by getMILPParams (optional) % % newConnected cell array with the reactions that could be % connected. This is not calulated if diff --git a/core/fitTasks.m b/core/fitTasks.m index 0a8cca96..9963f9e7 100755 --- a/core/fitTasks.m +++ b/core/fitTasks.m @@ -1,4 +1,4 @@ -function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params) +function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure) % fitTasks % Fills gaps in a model by including reactions from a reference model, % so that the resulting model can perform all the tasks in a task list @@ -17,7 +17,6 @@ % reactions (optional, default is -1 for all reactions) % taskStructure structure with the tasks, as from parseTaskList. If % this is supplied then inputFile is ignored (optional) -% params parameter structure as used by getMILPParams (optional) % % % Output: diff --git a/core/getMinNrFluxes.m b/core/getMinNrFluxes.m index 137e90b3..6c9ca5e2 100755 --- a/core/getMinNrFluxes.m +++ b/core/getMinNrFluxes.m @@ -8,7 +8,7 @@ % with the same number of elements as reactions in the model, % of a vector of indexes for the reactions that should be % minimized (optional, default model.rxns) -% params parameter structure as used by getMILPParams (optional) +% params *obsolete option* % scores vector of weights for the reactions. Negative scores % should not have flux. Positive scores are not possible in this % implementation, and they are changed to max(scores(scores<0)). diff --git a/core/makeSomething.m b/core/makeSomething.m index 4a5a753f..1b9b2ab8 100755 --- a/core/makeSomething.m +++ b/core/makeSomething.m @@ -21,7 +21,7 @@ % false) % allowExcretion allow for excretion of all other metabolites (optional, % default true) -% params parameter structure as used by getMILPParams (optional) +% params *obsolete option* % ignoreIntBounds true if internal bounds (including reversibility) % should be ignored. Exchange reactions are not affected. % This can be used to find unbalanced solutions which are diff --git a/core/randomSampling.m b/core/randomSampling.m index 0e41a308..9c48431c 100755 --- a/core/randomSampling.m +++ b/core/randomSampling.m @@ -42,13 +42,18 @@ % Output: % solutions matrix with the solutions % goodRxns double vector of indexes of those reactions -% that are not involved in loops and can be used -% as random objective functions +% that are not involved in loops or always carry +% zero flux and can be used as random objective +% functions % -% The solutions are generated by maximizing (with random weights) for a -% random set of three reactions. For reversible reactions it randomly +% Note: The solutions are generated by maximizing (with random weights) for +% a random set of three reactions. For reversible reactions it randomly % chooses between maximizing and minimizing. % +% If the model is a GECKO v3+ ecModel, then usage_prot reactions are not +% selected for sampling, instead focusing on sampling from the metabolic +% aspects that form the solution space. +% % Usage: solutions = randomSampling(model, nSamples, replaceBoundsWithInf,... % supressErrors, runParallel, goodRxns, minFlux) @@ -94,35 +99,39 @@ warning('The model objective function cannot reach a non-zero value. This might be intended, so randomSampling will continue, but this could indicate problems with your model') end -[~,hsSol]=solveLP(model); nRxns = numel(model.rxns); %Reactions which can be involved in loops should not be optimized for. %Check which reactions reach an arbitary high upper bound if nargin<6 || isempty(goodRxns) - goodRxns = true(numel(model.rxns),numel(model.rxns)); - goodRxns = num2cell(goodRxns,1); + revRxns = transpose(find(model.lb < 0)); + fwdRxns = transpose(find(model.ub > 0)); + rxnList = [fwdRxns,revRxns]; + objList = [ones(numel(fwdRxns),1);-ones(numel(revRxns),1)]; + testSol = zeros(numel(objList),1); PB = ProgressBar2(nRxns,'Prepare goodRxns not involved in loops','cli'); - parfor (i=1:nRxns) - testModel=setParam(model,'eq',model.rxns(i),1000); - sol=solveLP(testModel,0,[],hsSol); + + parfor i = 1:numel(objList) + testModel=setParam(model,'obj',rxnList(i),objList(i)); + sol=solveLP(testModel,0); if ~isempty(sol.f) - notGood = abs(sol.x)>999; - goodRxns{i}(notGood)=false; - else - %If the reaction is reversible, also check in that direction - if model.rev(i) - testModel=setParam(model,'eq',model.rxns(i),-1000); - sol=solveLP(testModel,0,[],hsSol); - if ~isempty(sol.f) - notGood = abs(sol.x)>999; - goodRxns{i}(notGood)=false; - end - end + testSol(i) = sol.x(rxnList(i)); end count(PB); end - goodRxns = cell2mat(goodRxns); - goodRxns = find(prod(goodRxns,2)); + testSol = testSol.*objList; + goodRxns = true(nRxns,1); + % Filter out reactions that can reach (-)1000 (= involved in loop) + goodRxns(fwdRxns(testSol(1:numel(fwdRxns)) > 999)) = false; + goodRxns(revRxns(testSol(numel(fwdRxns)+1:end) > 999)) = false; + testSol(revRxns) = testSol(revRxns) + testSol(numel(fwdRxns)+1:end); + % Filter out reactions that cannot carry flux + goodRxns(testSol(1:nRxns) == 0) = false; + % In ecModels do not sample from usage_prot reactions + if isfield(model,'ec') + usageRxns = startsWith(model.rxns,'usage_prot_'); + goodRxns(usageRxns) = false; + end + goodRxns = find(goodRxns); end %Reserve space for a solution matrix @@ -139,13 +148,16 @@ while lt(0,badSolutions) rxns = randsample(numel(goodRxns),nObjRxns); tmpModel.c = zeros(numel(tmpModel.rxns),1); + multipliers = randsample([-1 1],nObjRxns,true); - multipliers(tmpModel.rev(goodRxns(rxns))==0) = 1; + multipliers(tmpModel.ub(goodRxns(rxns))<0) = 1; + tmpModel.c(goodRxns(rxns)) = rand(nObjRxns,1).*multipliers; + if true(minFlux) - sol=solveLP(tmpModel,1,[],hsSol); + sol=solveLP(tmpModel,1); else - sol=solveLP(tmpModel,0,[],hsSol); + sol=solveLP(tmpModel,0); end if any(sol.x) && abs(sol.f)>10^-8 sols{i} = sol.x; diff --git a/core/setParam.m b/core/setParam.m index c341b914..e24ad2ca 100755 --- a/core/setParam.m +++ b/core/setParam.m @@ -123,4 +123,9 @@ end end end +if any(ismember(paramType,{'lb','ub','unc'})) + invalidBound = model.lb(indexes) > model.ub(indexes); + if any(invalidBound) + error(['Invalid set of bounds for reaction(s): ', strjoin(model.rxns(indexes(invalidBound)),', '), '.']) + end end diff --git a/doc/INIT/ftINIT.html b/doc/INIT/ftINIT.html index fed102da..3c05d6e0 100644 --- a/doc/INIT/ftINIT.html +++ b/doc/INIT/ftINIT.html @@ -82,9 +82,7 @@

DESCRIPTION ^SOURCE CODE ^% useScoresForTasks true if the calculated reaction scored should be 0055 % used as weights when fitting to tasks (optional, default 0056 % true) -0057 % paramsFT parameter structure as used by getMILPParams. This -0058 % is for the fitTasks step. For the INIT algorithm, -0059 % see params (optional, default []) -0060 % verbose if true, the MILP progression will be shown. -0061 % (optional, default false) -0062 % -0063 % model the resulting model structure -0064 % metProduction array that indicates which of the -0065 % metabolites in metabolomicsData that could be -0066 % produced. Note that this is before the -0067 % gap-filling process to enable defined tasks. To -0068 % see which metabolites that can be produced in -0069 % the final model, use canProduce. -0070 % -2: metabolite name not found in model -0071 % -1: metabolite found, but it could not be produced -0072 % 1: metabolite could be produced -0073 % addedRxnsForTasks cell array of the reactions which were added in -0074 % order to perform the tasks -0075 % deletedRxnsInINIT cell array of reactions deleted because they -0076 % could not carry flux (INIT requires a -0077 % functional input model) -0078 % fullMipRes The solver results from the last MILP step run -0079 % -0080 % This is the main function for automatic reconstruction of models based -0081 % on the ftINIT algorithm (). -0082 % -0083 % NOTE: Exchange metabolites should normally not be removed from the model -0084 % when using this approach, since checkTasks/fitTasks rely on putting specific -0085 % constraints for each task. The INIT algorithm will remove exchange metabolites -0086 % if any are present. Use importModel(file,false) to import a model with -0087 % exchange metabolites remaining. -0088 % -0089 % Usage: [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, ... -0090 % fullMipRes] = ... -0091 % ftINIT(prepData, tissue, celltype, hpaData, transcrData, ... -0092 % metabolomicsData, INITSteps, removeGenes, useScoresForTasks, ... -0093 % paramsFT); -0094 % -0095 -0096 -0097 if nargin < 5 -0098 transcrData = []; -0099 end -0100 if nargin < 6 -0101 metabolomicsData = []; -0102 end -0103 if nargin < 7 || isempty(INITSteps) -0104 INITSteps = getINITSteps([],'1+1'); -0105 end -0106 if nargin < 8 || isempty(removeGenes) -0107 removeGenes = true; -0108 end -0109 if nargin < 9 || isempty(useScoresForTasks) -0110 useScoresForTasks = true; -0111 end -0112 if nargin < 10 -0113 paramsFT = []; -0114 end -0115 -0116 if nargin < 11 -0117 verbose = false; -0118 end -0119 %Handle detected mets: -0120 %Previously, this was handled by giving a bonus for secreting those metabolites, -0121 %but that doesn't work since the metabolite secretion and uptake can be lost when -0122 %we merge linearly dependent reactions. -0123 %Instead, we need to figure out which reactions either produce or take up the mets. -0124 %We then give a bonus if any of them carry flux. -0125 %To simplify things, we focus on reactions that produce the metabolite (since there must be one such reaction). -0126 %It is still a bit complicated though. In this step, we focus on identifying -0127 %producer reactions. We further reason that the direction doesn't matter - -0128 %we can force one of these reactions in any direction - if it becomes a consumer, it will -0129 %automatically force another producer on as well (otherwise we'll have a net consumption). -0130 -0131 if (~isempty(metabolomicsData)) -0132 if length(unique(upper(metabolomicsData))) ~= length(metabolomicsData) -0133 dispEM('Metabolomics contains the same metabolite multiple times'); -0134 end -0135 metData = false(numel(metabolomicsData), length(prepData.minModel.rxns)); %one row per metabolite that is a boolean vector -0136 for i=1:numel(metabolomicsData) -0137 %Get the matching mets -0138 metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); -0139 prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers -0140 (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers -0141 %convert the production rxns from refModel to minModel -0142 prepData.groupIds; -0143 [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); -0144 grpIdsMerged = nan(length(prepData.minModel.rxns),1); -0145 grpIdsMerged(ia) = prepData.groupIds(ib); -0146 -0147 groupIdsPos = unique(prepData.groupIds(prodRxnsSel));%gets all group ids which includes a production rxn -0148 groupIdsPos = groupIdsPos(groupIdsPos ~= 0);%remove the 0 id, it means there is no group -0149 %the other option is that there is a direct match between the rxn id in minModel and refModel: -0150 posRxns = prepData.refModel.rxns(prodRxnsSel); -0151 directMatch = ismember(prepData.minModel.rxns, posRxns).'; -0152 -0153 metData(i,:) = ismember(grpIdsMerged, groupIdsPos).' | directMatch; -0154 end -0155 metData = sparse(metData); -0156 else -0157 metData = []; -0158 end -0159 -0160 % Get rxn scores and adapt them to the minimized model -0161 origRxnScores = scoreComplexModel(prepData.refModel,hpaData,transcrData,tissue,celltype); -0162 origRxnScores(origRxnScores > -0.1 & origRxnScores <= 0) = -0.1;%we don't want reaction scores that are exactly 0 (or close), this causes problems in the milp -0163 origRxnScores(origRxnScores < 0.1 & origRxnScores > 0) = 0.1; -0164 -0165 rxnsTurnedOn = false(length(prepData.minModel.rxns),1); -0166 fluxes = zeros(length(prepData.minModel.rxns),1); +0057 % paramsFT *obsolete option* +0058 % verbose if true, the MILP progression will be shown. +0059 % (optional, default false) +0060 % +0061 % model the resulting model structure +0062 % metProduction array that indicates which of the +0063 % metabolites in metabolomicsData that could be +0064 % produced. Note that this is before the +0065 % gap-filling process to enable defined tasks. To +0066 % see which metabolites that can be produced in +0067 % the final model, use canProduce. +0068 % -2: metabolite name not found in model +0069 % -1: metabolite found, but it could not be produced +0070 % 1: metabolite could be produced +0071 % addedRxnsForTasks cell array of the reactions which were added in +0072 % order to perform the tasks +0073 % deletedRxnsInINIT cell array of reactions deleted because they +0074 % could not carry flux (INIT requires a +0075 % functional input model) +0076 % fullMipRes The solver results from the last MILP step run +0077 % +0078 % This is the main function for automatic reconstruction of models based +0079 % on the ftINIT algorithm (). +0080 % +0081 % NOTE: Exchange metabolites should normally not be removed from the model +0082 % when using this approach, since checkTasks/fitTasks rely on putting specific +0083 % constraints for each task. The INIT algorithm will remove exchange metabolites +0084 % if any are present. Use importModel(file,false) to import a model with +0085 % exchange metabolites remaining. +0086 % +0087 % Usage: [model, metProduction, addedRxnsForTasks, deletedRxnsInINIT, ... +0088 % fullMipRes] = ... +0089 % ftINIT(prepData, tissue, celltype, hpaData, transcrData, ... +0090 % metabolomicsData, INITSteps, removeGenes, useScoresForTasks, ... +0091 % paramsFT); +0092 % +0093 +0094 +0095 if nargin < 5 +0096 transcrData = []; +0097 end +0098 if nargin < 6 +0099 metabolomicsData = []; +0100 end +0101 if nargin < 7 || isempty(INITSteps) +0102 INITSteps = getINITSteps([],'1+1'); +0103 end +0104 if nargin < 8 || isempty(removeGenes) +0105 removeGenes = true; +0106 end +0107 if nargin < 9 || isempty(useScoresForTasks) +0108 useScoresForTasks = true; +0109 end +0110 if nargin < 10 +0111 paramsFT = []; +0112 end +0113 +0114 if nargin < 11 +0115 verbose = false; +0116 end +0117 %Handle detected mets: +0118 %Previously, this was handled by giving a bonus for secreting those metabolites, +0119 %but that doesn't work since the metabolite secretion and uptake can be lost when +0120 %we merge linearly dependent reactions. +0121 %Instead, we need to figure out which reactions either produce or take up the mets. +0122 %We then give a bonus if any of them carry flux. +0123 %To simplify things, we focus on reactions that produce the metabolite (since there must be one such reaction). +0124 %It is still a bit complicated though. In this step, we focus on identifying +0125 %producer reactions. We further reason that the direction doesn't matter - +0126 %we can force one of these reactions in any direction - if it becomes a consumer, it will +0127 %automatically force another producer on as well (otherwise we'll have a net consumption). +0128 +0129 if (~isempty(metabolomicsData)) +0130 if length(unique(upper(metabolomicsData))) ~= length(metabolomicsData) +0131 dispEM('Metabolomics contains the same metabolite multiple times'); +0132 end +0133 metData = false(numel(metabolomicsData), length(prepData.minModel.rxns)); %one row per metabolite that is a boolean vector +0134 for i=1:numel(metabolomicsData) +0135 %Get the matching mets +0136 metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); +0137 prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers +0138 (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers +0139 %convert the production rxns from refModel to minModel +0140 prepData.groupIds; +0141 [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); +0142 grpIdsMerged = nan(length(prepData.minModel.rxns),1); +0143 grpIdsMerged(ia) = prepData.groupIds(ib); +0144 +0145 groupIdsPos = unique(prepData.groupIds(prodRxnsSel));%gets all group ids which includes a production rxn +0146 groupIdsPos = groupIdsPos(groupIdsPos ~= 0);%remove the 0 id, it means there is no group +0147 %the other option is that there is a direct match between the rxn id in minModel and refModel: +0148 posRxns = prepData.refModel.rxns(prodRxnsSel); +0149 directMatch = ismember(prepData.minModel.rxns, posRxns).'; +0150 +0151 metData(i,:) = ismember(grpIdsMerged, groupIdsPos).' | directMatch; +0152 end +0153 metData = sparse(metData); +0154 else +0155 metData = []; +0156 end +0157 +0158 % Get rxn scores and adapt them to the minimized model +0159 origRxnScores = scoreComplexModel(prepData.refModel,hpaData,transcrData,tissue,celltype); +0160 origRxnScores(origRxnScores > -0.1 & origRxnScores <= 0) = -0.1;%we don't want reaction scores that are exactly 0 (or close), this causes problems in the milp +0161 origRxnScores(origRxnScores < 0.1 & origRxnScores > 0) = 0.1; +0162 +0163 rxnsTurnedOn = false(length(prepData.minModel.rxns),1); +0164 fluxes = zeros(length(prepData.minModel.rxns),1); +0165 +0166 rxnsToIgnoreLastStep = [1;1;1;1;1;1;1;1]; 0167 -0168 rxnsToIgnoreLastStep = [1;1;1;1;1;1;1;1]; -0169 -0170 %We assume that all essential rxns are irrev - this is taken care of in -0171 %prepINITModel. We then use an initial flux "from last run" of 0.1 for all -0172 %reactions. This is used for knowing what flux should be forced through an -0173 %essential rxn. -0174 fluxes = ones(length(prepData.minModel.rxns), 1).*0.1; -0175 -0176 for initStep = 1:length(INITSteps) -0177 disp(['ftINIT: Running step ' num2str(initStep)]) -0178 stp = INITSteps{initStep}; -0179 -0180 if any ((rxnsToIgnoreLastStep - stp.RxnsToIgnoreMask) < 0) -0181 dispEM('RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.'); -0182 end -0183 rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask; +0168 %We assume that all essential rxns are irrev - this is taken care of in +0169 %prepINITModel. We then use an initial flux "from last run" of 0.1 for all +0170 %reactions. This is used for knowing what flux should be forced through an +0171 %essential rxn. +0172 fluxes = ones(length(prepData.minModel.rxns), 1).*0.1; +0173 +0174 for initStep = 1:length(INITSteps) +0175 disp(['ftINIT: Running step ' num2str(initStep)]) +0176 stp = INITSteps{initStep}; +0177 +0178 if any ((rxnsToIgnoreLastStep - stp.RxnsToIgnoreMask) < 0) +0179 dispEM('RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.'); +0180 end +0181 rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask; +0182 +0183 mm = prepData.minModel; 0184 -0185 mm = prepData.minModel; -0186 -0187 if (~isempty(stp.MetsToIgnore)) -0188 if (~isempty(stp.MetsToIgnore.simpleMets)) -0189 %Here, we remove simple metabolites that will not really affect the milp but -0190 %are very common in the S matrix. For example H2O, H+, etc. -0191 %It is also possible to leave compartments untouched, for example the i compartment in the mitochondria (for H+). -0192 metsToRem = ismember(mm.metNames,stp.MetsToIgnore.simpleMets.mets); -0193 compsToKeep = find(ismember(mm.comps, stp.MetsToIgnore.simpleMets.compsToKeep)); -0194 metsToRem = metsToRem & ~ismember(mm.metComps, compsToKeep); -0195 mm.S(metsToRem,:) = 0; -0196 end -0197 end -0198 -0199 %Set up the reaction scores and essential rxns -0200 rxnsToIgnore = getRxnsFromPattern(stp.RxnsToIgnoreMask, prepData); -0201 rxnScores = groupRxnScores(prepData.minModel, origRxnScores, prepData.refModel.rxns, prepData.groupIds, rxnsToIgnore); -0202 -0203 essentialRxns = prepData.essentialRxns; -0204 toRev = false(numel(mm.rxns),1); -0205 %Handle the results from previous steps ('ignore', 'exclude', 'essential') -0206 if strcmp(stp.HowToUsePrevResults, 'exclude') -0207 rxnScores(rxnsTurnedOn) = 0; %This is not used anymore in any step setup. -0208 elseif strcmp(stp.HowToUsePrevResults, 'essential') -0209 %Make all reversible reactions turned on in previous steps reversible -0210 %in the direction that they were previously carrying flux -0211 -0212 %first reverse the reactions that need to be reversed -0213 rev = mm.rev == 1; -0214 toRev = rxnsTurnedOn & rev & fluxes < 0; -0215 mm = reverseRxns(mm, mm.rxns(toRev)); -0216 -0217 %Then make them irreversible -0218 mm.rev(rxnsTurnedOn) = 0; -0219 mm.lb(rxnsTurnedOn) = 0; -0220 -0221 essentialRxns = unique([prepData.essentialRxns;mm.rxns(rxnsTurnedOn)]); -0222 end -0223 -0224 -0225 mipGap = 1; -0226 first = true; -0227 success = false; -0228 fullMipRes = []; -0229 for rn = 1:length(stp.MILPParams) -0230 params = stp.MILPParams{rn}; -0231 if ~isfield(params, 'MIPGap') -0232 params.MIPGap = 0.0004; -0233 end -0234 -0235 if ~isfield(params, 'TimeLimit') -0236 params.TimeLimit = 5000; -0237 end -0238 -0239 if ~first -0240 %There is sometimes a problem with that the objective function becomes close to zero, -0241 %which leads to that a small percentage of that (which is the MIPGap sent in) is very small -0242 %and the MILP hence takes a lot of time to finish. We also therefore use an absolute MIP gap, -0243 %converted to a percentage using the last value of the objective function. -0244 params.MIPGap = min(max(params.MIPGap, stp.AbsMIPGaps{rn}/abs(lastObjVal)),1); -0245 params.seed = 1234;%use another seed, may work better -0246 -0247 if mipGap <= params.MIPGap -0248 success = true; -0249 break; %we're done - this will not happen the first time -0250 else -0251 disp(['MipGap too high, trying with a different run. MipGap = ' num2str(mipGap) ' New MipGap Limit = ' num2str(params.MIPGap)]) -0252 end -0253 end +0185 if (~isempty(stp.MetsToIgnore)) +0186 if (~isempty(stp.MetsToIgnore.simpleMets)) +0187 %Here, we remove simple metabolites that will not really affect the milp but +0188 %are very common in the S matrix. For example H2O, H+, etc. +0189 %It is also possible to leave compartments untouched, for example the i compartment in the mitochondria (for H+). +0190 metsToRem = ismember(mm.metNames,stp.MetsToIgnore.simpleMets.mets); +0191 compsToKeep = find(ismember(mm.comps, stp.MetsToIgnore.simpleMets.compsToKeep)); +0192 metsToRem = metsToRem & ~ismember(mm.metComps, compsToKeep); +0193 mm.S(metsToRem,:) = 0; +0194 end +0195 end +0196 +0197 %Set up the reaction scores and essential rxns +0198 rxnsToIgnore = getRxnsFromPattern(stp.RxnsToIgnoreMask, prepData); +0199 rxnScores = groupRxnScores(prepData.minModel, origRxnScores, prepData.refModel.rxns, prepData.groupIds, rxnsToIgnore); +0200 +0201 essentialRxns = prepData.essentialRxns; +0202 toRev = false(numel(mm.rxns),1); +0203 %Handle the results from previous steps ('ignore', 'exclude', 'essential') +0204 if strcmp(stp.HowToUsePrevResults, 'exclude') +0205 rxnScores(rxnsTurnedOn) = 0; %This is not used anymore in any step setup. +0206 elseif strcmp(stp.HowToUsePrevResults, 'essential') +0207 %Make all reversible reactions turned on in previous steps reversible +0208 %in the direction that they were previously carrying flux +0209 +0210 %first reverse the reactions that need to be reversed +0211 rev = mm.rev == 1; +0212 toRev = rxnsTurnedOn & rev & fluxes < 0; +0213 mm = reverseRxns(mm, mm.rxns(toRev)); +0214 +0215 %Then make them irreversible +0216 mm.rev(rxnsTurnedOn) = 0; +0217 mm.lb(rxnsTurnedOn) = 0; +0218 +0219 essentialRxns = unique([prepData.essentialRxns;mm.rxns(rxnsTurnedOn)]); +0220 end +0221 +0222 +0223 mipGap = 1; +0224 first = true; +0225 success = false; +0226 fullMipRes = []; +0227 for rn = 1:length(stp.MILPParams) +0228 params = stp.MILPParams{rn}; +0229 if ~isfield(params, 'MIPGap') +0230 params.MIPGap = 0.0004; +0231 end +0232 +0233 if ~isfield(params, 'TimeLimit') +0234 params.TimeLimit = 5000; +0235 end +0236 +0237 if ~first +0238 %There is sometimes a problem with that the objective function becomes close to zero, +0239 %which leads to that a small percentage of that (which is the MIPGap sent in) is very small +0240 %and the MILP hence takes a lot of time to finish. We also therefore use an absolute MIP gap, +0241 %converted to a percentage using the last value of the objective function. +0242 params.MIPGap = min(max(params.MIPGap, stp.AbsMIPGaps{rn}/abs(lastObjVal)),1); +0243 params.seed = 1234;%use another seed, may work better +0244 +0245 if mipGap <= params.MIPGap +0246 success = true; +0247 break; %we're done - this will not happen the first time +0248 else +0249 disp(['MipGap too high, trying with a different run. MipGap = ' num2str(mipGap) ' New MipGap Limit = ' num2str(params.MIPGap)]) +0250 end +0251 end +0252 +0253 first = false; 0254 -0255 first = false; -0256 -0257 %now run the MILP -0258 try -0259 %The prodweight for metabolomics is currently set to 5 - 0.5 was default in the old version, which I deemed very small? -0260 %There could be a need to specify this somewhere in the call at some point. -0261 %This value has not been evaluated, but is assumed in the test cases - if changed, update the test case -0262 startVals = []; -0263 if ~isempty(fullMipRes) -0264 startVals = fullMipRes.full; -0265 end -0266 [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose); -0267 %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed -0268 fluxes1(toRev) = -fluxes1(toRev); -0269 -0270 mipGap = fullMipRes.mipgap; -0271 lastObjVal = fullMipRes.obj; -0272 catch e -0273 mipGap = Inf; -0274 lastObjVal = Inf; %we need to set something here, Inf leads to that this doesn't come into play -0275 end -0276 -0277 success = mipGap <= params.MIPGap; -0278 end -0279 -0280 if ~success -0281 dispEM(['Failed to find good enough solution within the time frame. MIPGap: ' num2str(mipGap)]); -0282 end -0283 -0284 %save the reactions turned on and their fluxes for the next step -0285 rxnsTurnedOn = rxnsTurnedOn | rxnsTurnedOn1.'; -0286 %The fluxes are a bit tricky - what if they change direction between the steps? -0287 %The fluxes are used to determine the direction in which reactions are forced on -0288 %(to simplify the problem it is good if they are unidirectional). -0289 %We use the following strategy: -0290 %1. Use the fluxes from the most recent step. -0291 %2. If any flux is very low there (i.e. basically zero), use the flux from the previous steps -0292 %This could in theory cause problems, but seems to work well practically -0293 fluxesOld = fluxes; -0294 fluxes = fluxes1; -0295 %make sure that all reactions that are on actually has a flux - otherwise -0296 %things could go bad, since the flux will be set to essential in a random direction -0297 %This sometimes happens for rxns with negative score - let's just accept that. -0298 %if (sum(abs(fluxes1) < 10^-7 & rxnsTurnedOn)) -0299 % dispEM('There are rxns turned on without flux - this might cause problems'); -0300 %end -0301 %fluxes(abs(fluxes1) < 10^-7) = fluxesOld(abs(fluxes1) < 10^-9); -0302 end -0303 -0304 -0305 %get the essential rxns -0306 essential = ismember(prepData.minModel.rxns,prepData.essentialRxns); -0307 %So, we only add reactions where the linearly merged scores are zero for all linearly dependent reactions -0308 % (this cannot happen by chance, taken care of in the function groupRxnScores) -0309 rxnsToIgn = rxnScores == 0; -0310 deletedRxnsInINITSel = ~(rxnsTurnedOn | rxnsToIgn | essential); -0311 deletedRxnsInINIT = prepData.minModel.rxns(deletedRxnsInINITSel); -0312 -0313 %Here we need to figure out which original reactions (before the linear merge) -0314 %that were removed. These are all reactions with the same group ids as the removed reactions -0315 groupIdsRemoved = prepData.groupIds(ismember(prepData.refModel.rxns, deletedRxnsInINIT)); %can improve this slightly, use sel above -0316 groupIdsRemoved = groupIdsRemoved(groupIdsRemoved ~= 0);%zero means that the reaction was not grouped, all with zeros are not a group! -0317 rxnsToRem = union(prepData.refModel.rxns(ismember(prepData.groupIds,groupIdsRemoved)), deletedRxnsInINIT);%make a union here to include the ungrouped (unmerged) as well +0255 %now run the MILP +0256 try +0257 %The prodweight for metabolomics is currently set to 5 - 0.5 was default in the old version, which I deemed very small? +0258 %There could be a need to specify this somewhere in the call at some point. +0259 %This value has not been evaluated, but is assumed in the test cases - if changed, update the test case +0260 startVals = []; +0261 if ~isempty(fullMipRes) +0262 startVals = fullMipRes.full; +0263 end +0264 [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose); +0265 %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed +0266 fluxes1(toRev) = -fluxes1(toRev); +0267 +0268 mipGap = fullMipRes.mipgap; +0269 lastObjVal = fullMipRes.obj; +0270 catch e +0271 mipGap = Inf; +0272 lastObjVal = Inf; %we need to set something here, Inf leads to that this doesn't come into play +0273 end +0274 +0275 success = mipGap <= params.MIPGap; +0276 end +0277 +0278 if ~success +0279 dispEM(['Failed to find good enough solution within the time frame. MIPGap: ' num2str(mipGap)]); +0280 end +0281 +0282 %save the reactions turned on and their fluxes for the next step +0283 rxnsTurnedOn = rxnsTurnedOn | rxnsTurnedOn1.'; +0284 %The fluxes are a bit tricky - what if they change direction between the steps? +0285 %The fluxes are used to determine the direction in which reactions are forced on +0286 %(to simplify the problem it is good if they are unidirectional). +0287 %We use the following strategy: +0288 %1. Use the fluxes from the most recent step. +0289 %2. If any flux is very low there (i.e. basically zero), use the flux from the previous steps +0290 %This could in theory cause problems, but seems to work well practically +0291 fluxesOld = fluxes; +0292 fluxes = fluxes1; +0293 %make sure that all reactions that are on actually has a flux - otherwise +0294 %things could go bad, since the flux will be set to essential in a random direction +0295 %This sometimes happens for rxns with negative score - let's just accept that. +0296 %if (sum(abs(fluxes1) < 10^-7 & rxnsTurnedOn)) +0297 % dispEM('There are rxns turned on without flux - this might cause problems'); +0298 %end +0299 %fluxes(abs(fluxes1) < 10^-7) = fluxesOld(abs(fluxes1) < 10^-9); +0300 end +0301 +0302 +0303 %get the essential rxns +0304 essential = ismember(prepData.minModel.rxns,prepData.essentialRxns); +0305 %So, we only add reactions where the linearly merged scores are zero for all linearly dependent reactions +0306 % (this cannot happen by chance, taken care of in the function groupRxnScores) +0307 rxnsToIgn = rxnScores == 0; +0308 deletedRxnsInINITSel = ~(rxnsTurnedOn | rxnsToIgn | essential); +0309 deletedRxnsInINIT = prepData.minModel.rxns(deletedRxnsInINITSel); +0310 +0311 %Here we need to figure out which original reactions (before the linear merge) +0312 %that were removed. These are all reactions with the same group ids as the removed reactions +0313 groupIdsRemoved = prepData.groupIds(ismember(prepData.refModel.rxns, deletedRxnsInINIT)); %can improve this slightly, use sel above +0314 groupIdsRemoved = groupIdsRemoved(groupIdsRemoved ~= 0);%zero means that the reaction was not grouped, all with zeros are not a group! +0315 rxnsToRem = union(prepData.refModel.rxns(ismember(prepData.groupIds,groupIdsRemoved)), deletedRxnsInINIT);%make a union here to include the ungrouped (unmerged) as well +0316 +0317 initModel = removeReactions(prepData.refModel,rxnsToRem,false,true); 0318 -0319 initModel = removeReactions(prepData.refModel,rxnsToRem,false,true); -0320 -0321 % remove metabolites separately to avoid removing those needed for tasks -0322 unusedMets = initModel.mets(all(initModel.S == 0,2)); -0323 initModel = removeMets(initModel, setdiff(unusedMets, prepData.essentialMetsForTasks)); -0324 -0325 %if printReport == true -0326 % printScores(initModel,'INIT model statistics',hpaData,transcrData,tissue,celltype); -0327 % printScores(removeReactions(cModel,setdiff(cModel.rxns,rxnsToRem),true,true),'Reactions deleted by INIT',hpaData,transcrData,tissue,celltype); -0328 %end -0329 -0330 %The full model has exchange reactions in it. ftINITFillGapsForAllTasks calls -0331 %ftINITFillGaps, which automatically removes exchange metabolites (because it -0332 %assumes that the reactions are constrained when appropriate). In this case the -0333 %uptakes/outputs are retrieved from the task sheet instead. To prevent -0334 %exchange reactions being used to fill gaps, they are deleted from the -0335 %reference model here. -0336 initModel.id = 'INITModel'; -0337 -0338 %If gaps in the model should be filled using a task list -0339 if ~isempty(prepData.taskStruct) -0340 %Remove exchange reactions and reactions already included in the INIT -0341 %model -0342 %We changed strategy and instead include all rxns except the exchange rxns in the ref model -0343 %But we do keep the exchange rxns that are essential. -0344 %Let's test to remove all, that should work -0345 -0346 %At this stage the model is fully connected and most of the genes with -0347 %good scores should have been included. The final gap-filling should -0348 %take the scores of the genes into account, so that "rather bad" -0349 %reactions are preferred to "very bad" reactions. However, reactions -0350 %with positive scores will be included even if they are not connected -0351 %in the current formulation. Therefore, such reactions will have to be -0352 %assigned a small negative score instead. -0353 exchRxns = getExchangeRxns(prepData.refModel); -0354 refModelNoExc = removeReactions(prepData.refModelWithBM,exchRxns,false,true); -0355 exchRxns = getExchangeRxns(initModel); -0356 initModelNoExc = removeReactions(closeModel(initModel),exchRxns,false,true); -0357 -0358 if useScoresForTasks == true -0359 %map the rxn scores to the model without exchange rxns -0360 [~,ia,ib] = intersect(refModelNoExc.rxns,prepData.refModel.rxns); -0361 rxnScores2nd = NaN(length(refModelNoExc.rxns),1); -0362 rxnScores2nd(ia) = origRxnScores(ib); -0363 %all(rxnScores2nd == refRxnScores);%should be the same, ok! -0364 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,min(rxnScores2nd,-0.1),prepData.taskStruct,paramsFT,verbose); -0365 else -0366 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,[],prepData.taskStruct,paramsFT,verbose); -0367 end -0368 %if printReport == true -0369 % printScores(outModel,'Functional model statistics',hpaData,transcrData,tissue,celltype); -0370 % printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,transcrData,tissue,celltype); -0371 %end -0372 -0373 addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2)); -0374 else -0375 outModel = initModel; -0376 addedRxnMat = []; -0377 addedRxnsForTasks = {}; -0378 end -0379 -0380 % The model can now perform all the tasks defined in the task list. -0381 model = outModel; -0382 -0383 -0384 % At this stage the model will contain some exchange reactions but probably -0385 % not all (and maybe zero). This can be inconvenient, so all exchange -0386 % reactions from the reference model are added, except for those which -0387 % involve metabolites that are not in the model. -0388 -0389 %Start from the original model, and just remove the reactions that are no longer there (and keep exchange rxns). The model we got out -0390 %from the problem is not complete, it doesn't have GRPs etc. -0391 %The logic below is a bit complicated. We identify the reactions that should be removed from the full model as -0392 %reactions that have been removed in the init model except the ones that were added back. In addition, we make -0393 %sure that no exchange rxns are removed - they can be removed in the init model if they were linearly merged with other -0394 %reactions that were decided to be removed from the model. We want to keep all exchange rxns to make sure the tasks can -0395 %be performed also without manipulating the b vector in the model (which is what is done in the gap-filling). -0396 exchRxns = getExchangeRxns(prepData.refModel); -0397 deletedRxnsInINIT = setdiff(prepData.refModel.rxns,union(union(initModel.rxns, addedRxnsForTasks), exchRxns)); -0398 outModel = removeReactions(prepData.refModel, deletedRxnsInINIT, true); %we skip removing the genes for now, I'm not sure it is desirable -0399 -0400 % If requested, attempt to remove negative-score genes from the model, -0401 % depending on their role (isozyme or complex subunit) in each grRule. -0402 % See the "removeLowScoreGenes" function more more details, and to adjust -0403 % any default parameters therein. -0404 if ( removeGenes ) -0405 [~, geneScores] = scoreComplexModel(outModel,hpaData,transcrData,tissue,celltype); -0406 outModel = removeLowScoreGenes(outModel,geneScores); -0407 end -0408 +0319 % remove metabolites separately to avoid removing those needed for tasks +0320 unusedMets = initModel.mets(all(initModel.S == 0,2)); +0321 initModel = removeMets(initModel, setdiff(unusedMets, prepData.essentialMetsForTasks)); +0322 +0323 %if printReport == true +0324 % printScores(initModel,'INIT model statistics',hpaData,transcrData,tissue,celltype); +0325 % printScores(removeReactions(cModel,setdiff(cModel.rxns,rxnsToRem),true,true),'Reactions deleted by INIT',hpaData,transcrData,tissue,celltype); +0326 %end +0327 +0328 %The full model has exchange reactions in it. ftINITFillGapsForAllTasks calls +0329 %ftINITFillGaps, which automatically removes exchange metabolites (because it +0330 %assumes that the reactions are constrained when appropriate). In this case the +0331 %uptakes/outputs are retrieved from the task sheet instead. To prevent +0332 %exchange reactions being used to fill gaps, they are deleted from the +0333 %reference model here. +0334 initModel.id = 'INITModel'; +0335 +0336 %If gaps in the model should be filled using a task list +0337 if ~isempty(prepData.taskStruct) +0338 %Remove exchange reactions and reactions already included in the INIT +0339 %model +0340 %We changed strategy and instead include all rxns except the exchange rxns in the ref model +0341 %But we do keep the exchange rxns that are essential. +0342 %Let's test to remove all, that should work +0343 +0344 %At this stage the model is fully connected and most of the genes with +0345 %good scores should have been included. The final gap-filling should +0346 %take the scores of the genes into account, so that "rather bad" +0347 %reactions are preferred to "very bad" reactions. However, reactions +0348 %with positive scores will be included even if they are not connected +0349 %in the current formulation. Therefore, such reactions will have to be +0350 %assigned a small negative score instead. +0351 exchRxns = getExchangeRxns(prepData.refModel); +0352 refModelNoExc = removeReactions(prepData.refModelWithBM,exchRxns,false,true); +0353 exchRxns = getExchangeRxns(initModel); +0354 initModelNoExc = removeReactions(closeModel(initModel),exchRxns,false,true); +0355 +0356 if useScoresForTasks == true +0357 %map the rxn scores to the model without exchange rxns +0358 [~,ia,ib] = intersect(refModelNoExc.rxns,prepData.refModel.rxns); +0359 rxnScores2nd = NaN(length(refModelNoExc.rxns),1); +0360 rxnScores2nd(ia) = origRxnScores(ib); +0361 %all(rxnScores2nd == refRxnScores);%should be the same, ok! +0362 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,min(rxnScores2nd,-0.1),prepData.taskStruct,paramsFT,verbose); +0363 else +0364 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,[],prepData.taskStruct,paramsFT,verbose); +0365 end +0366 %if printReport == true +0367 % printScores(outModel,'Functional model statistics',hpaData,transcrData,tissue,celltype); +0368 % printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,transcrData,tissue,celltype); +0369 %end +0370 +0371 addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2)); +0372 else +0373 outModel = initModel; +0374 addedRxnMat = []; +0375 addedRxnsForTasks = {}; +0376 end +0377 +0378 % The model can now perform all the tasks defined in the task list. +0379 model = outModel; +0380 +0381 +0382 % At this stage the model will contain some exchange reactions but probably +0383 % not all (and maybe zero). This can be inconvenient, so all exchange +0384 % reactions from the reference model are added, except for those which +0385 % involve metabolites that are not in the model. +0386 +0387 %Start from the original model, and just remove the reactions that are no longer there (and keep exchange rxns). The model we got out +0388 %from the problem is not complete, it doesn't have GRPs etc. +0389 %The logic below is a bit complicated. We identify the reactions that should be removed from the full model as +0390 %reactions that have been removed in the init model except the ones that were added back. In addition, we make +0391 %sure that no exchange rxns are removed - they can be removed in the init model if they were linearly merged with other +0392 %reactions that were decided to be removed from the model. We want to keep all exchange rxns to make sure the tasks can +0393 %be performed also without manipulating the b vector in the model (which is what is done in the gap-filling). +0394 exchRxns = getExchangeRxns(prepData.refModel); +0395 deletedRxnsInINIT = setdiff(prepData.refModel.rxns,union(union(initModel.rxns, addedRxnsForTasks), exchRxns)); +0396 outModel = removeReactions(prepData.refModel, deletedRxnsInINIT, true); %we skip removing the genes for now, I'm not sure it is desirable +0397 +0398 % If requested, attempt to remove negative-score genes from the model, +0399 % depending on their role (isozyme or complex subunit) in each grRule. +0400 % See the "removeLowScoreGenes" function more more details, and to adjust +0401 % any default parameters therein. +0402 if ( removeGenes ) +0403 [~, geneScores] = scoreComplexModel(outModel,hpaData,transcrData,tissue,celltype); +0404 outModel = removeLowScoreGenes(outModel,geneScores); +0405 end +0406 +0407 +0408 model = outModel; 0409 -0410 model = outModel; +0410 end 0411 -0412 end -0413 -0414 %This is for printing a summary of a model -0415 function [rxnS, geneS] = printScores(model,name,hpaData,transcrData,tissue,celltype) -0416 [a, b] = scoreComplexModel(model,hpaData,transcrData,tissue,celltype); -0417 rxnS = mean(a); -0418 geneS = mean(b,'omitnan'); -0419 fprintf([name ':\n']); -0420 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']); -0421 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']); -0422 fprintf(['\tMean gene score: ' num2str(geneS) '\n']); -0423 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']); -0424 end -0425 -0426 function rxnsToIgnore = getRxnsFromPattern(rxnsToIgnorePattern, prepData) -0427 rxnsToIgnore = false(length(prepData.toIgnoreExch),1); -0428 if rxnsToIgnorePattern(1) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreExch; end; -0429 if rxnsToIgnorePattern(2) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreImportRxns; end; -0430 if rxnsToIgnorePattern(3) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSimpleTransp; end; -0431 if rxnsToIgnorePattern(4) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAdvTransp; end; -0432 if rxnsToIgnorePattern(5) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSpont; end; -0433 if rxnsToIgnorePattern(6) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreS; end; -0434 if rxnsToIgnorePattern(7) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreCustomRxns; end; -0435 if rxnsToIgnorePattern(8) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAllWithoutGPRs; end; -0436 end -0437 +0412 %This is for printing a summary of a model +0413 function [rxnS, geneS] = printScores(model,name,hpaData,transcrData,tissue,celltype) +0414 [a, b] = scoreComplexModel(model,hpaData,transcrData,tissue,celltype); +0415 rxnS = mean(a); +0416 geneS = mean(b,'omitnan'); +0417 fprintf([name ':\n']); +0418 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']); +0419 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']); +0420 fprintf(['\tMean gene score: ' num2str(geneS) '\n']); +0421 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']); +0422 end +0423 +0424 function rxnsToIgnore = getRxnsFromPattern(rxnsToIgnorePattern, prepData) +0425 rxnsToIgnore = false(length(prepData.toIgnoreExch),1); +0426 if rxnsToIgnorePattern(1) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreExch; end; +0427 if rxnsToIgnorePattern(2) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreImportRxns; end; +0428 if rxnsToIgnorePattern(3) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSimpleTransp; end; +0429 if rxnsToIgnorePattern(4) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAdvTransp; end; +0430 if rxnsToIgnorePattern(5) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSpont; end; +0431 if rxnsToIgnorePattern(6) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreS; end; +0432 if rxnsToIgnorePattern(7) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreCustomRxns; end; +0433 if rxnsToIgnorePattern(8) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAllWithoutGPRs; end; +0434 end +0435
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/INIT/ftINITFillGaps.html b/doc/INIT/ftINITFillGaps.html index 6cefba5d..6a6661e7 100644 --- a/doc/INIT/ftINITFillGaps.html +++ b/doc/INIT/ftINITFillGaps.html @@ -51,7 +51,7 @@

DESCRIPTION ^SOURCE CODE ^% reference tModel. 0024 % The solver will try to maximize the sum of the 0025 % scores for the included reactions -0026 % params parameter structure as used by getMILPParams +0026 % params *obsolete option* 0027 % verbose if true, the MILP progression will be shown. 0028 % 0029 % addedRxns the rxns added diff --git a/doc/INIT/ftINITFillGapsForAllTasks.html b/doc/INIT/ftINITFillGapsForAllTasks.html index a234f544..ae7fdaf2 100644 --- a/doc/INIT/ftINITFillGapsForAllTasks.html +++ b/doc/INIT/ftINITFillGapsForAllTasks.html @@ -46,7 +46,7 @@

DESCRIPTION ^SOURCE CODE ^% reactions (optional, default is -1 for all reactions) 0019 % taskStructure structure with the tasks, as from parseTaskList. If 0020 % this is supplied then inputFile is ignored -0021 % params parameter structure as used by getMILPParams +0021 % params *obsolete option* 0022 % verbose if true, the MILP progression will be shown. 0023 % 0024 % diff --git a/doc/INIT/ftINITFillGapsMILP.html b/doc/INIT/ftINITFillGapsMILP.html index 6dcefd97..f2294249 100644 --- a/doc/INIT/ftINITFillGapsMILP.html +++ b/doc/INIT/ftINITFillGapsMILP.html @@ -40,7 +40,7 @@

DESCRIPTION ^SOURCE CODE ^% with the same number of elements as reactions in the model, 0013 % of a vector of indexes for the reactions that should be 0014 % minimized (optional, default model.rxns) -0015 % params parameter structure as used by getMILPParams (optional) +0015 % params *obsolete option* 0016 % scores vector of weights for the reactions. Negative scores 0017 % should not have flux. Positive scores are not possible in this 0018 % implementation, and they are changed to max(scores(scores<0)). diff --git a/doc/INIT/getINITModel.html b/doc/INIT/getINITModel.html index a1079434..739882c2 100644 --- a/doc/INIT/getINITModel.html +++ b/doc/INIT/getINITModel.html @@ -24,7 +24,7 @@

PURPOSE ^getINITModel_legacy

SYNOPSIS ^

-
function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure, params, paramsFT)
+
function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure)

DESCRIPTION ^

 getINITModel_legacy
@@ -32,7 +32,7 @@ 
DESCRIPTION ^DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

 This function calls:
 
+
  • runINIT	runINIT
    • 
 This function is called by:
 
      
 
    
@@ -151,13 +144,13 @@ 
    SUBFUNCTIONS ^function [rxnS, geneS]=printScores(model,name,hpaData,arrayData,tissue,celltype)
 
 
    SOURCE CODE ^
    
-
    0001 function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure, params, paramsFT)
+
    0001 function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure)
 0002 % getINITModel_legacy
 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 original
 0005 %   implementation of tINIT, which is replaced by ftINIT.
 0006 %
-0007 %   Input:
+0007 % Input:
 0008 %   refModel            a model structure. The model should be in the
 0009 %                       closed form (no exchange reactions open). Import
 0010 %                       using import(filename,false). If the model is not
@@ -206,418 +199,411 @@ 
    SOURCE CODE ^%                       as an alternative way to define tasks when Excel
 0054 %                       sheets are not suitable. Overrides taskFile (optional,
 0055 %                       default [])
-0056 %   params              parameter structure as used by getMILPParams. This is
-0057 %                       for the INIT algorithm. For the the MILP problems
-0058 %                       solved to fit tasks, see paramsFT (optional, default [])
-0059 %   paramsFT            parameter structure as used by getMILPParams. This is
-0060 %                       for the fitTasks step. For the INIT algorithm, see
-0061 %                       params (optional, default [])
-0062 %
-0063 %
-0064 %   Output:
-0065 %   model                   the resulting model structure
-0066 %   metProduction           array that indicates which of the
-0067 %                           metabolites in metabolomicsData that could be
-0068 %                           produced. Note that this is before the
-0069 %                           gap-filling process to enable defined tasks. To
-0070 %                           see which metabolites that can be produced in
-0071 %                           the final model, use canProduce.
-0072 %                           -2: metabolite name not found in model
-0073 %                           -1: metabolite found, but it could not be produced
-0074 %                           1: metabolite could be produced
-0075 %   essentialRxnsForTasks   cell array of the reactions which were
-0076 %                           essential to perform the tasks
-0077 %   addedRxnsForTasks       cell array of the reactions which were added in
-0078 %                           order to perform the tasks
-0079 %   deletedDeadEndRxns      cell array of reactions deleted because they
-0080 %                           could not carry flux (INIT requires a
-0081 %                           functional input model)
-0082 %   deletedRxnsInINIT       cell array of the reactions which were deleted by
-0083 %                           the INIT algorithm
-0084 %   taskReport              structure with the results for each task
-0085 %       id                  cell array with the id of the task
-0086 %       description         cell array with the description of the task
-0087 %       ok                  boolean array with true if the task was successful
-0088 %       essential           cell array with cell arrays of essential
-0089 %                           reactions for the task
-0090 %       gapfill             cell array of cell arrays of reactions included
-0091 %                           in the gap-filling for the task
-0092 %
-0093 %   This is the main function for automatic reconstruction of models based
-0094 %   on the (t)INIT algorithm (PLoS Comput Biol. 2012;8(5):e1002518,
-0095 %   Mol Syst Biol. 2014;10:721). Not all settings are possible using this
-0096 %   function, and you may want to call the functions scoreModel, runINIT
-0097 %   and fitTasks individually instead.
-0098 %
-0099 %   NOTE: Exchange metabolites should normally not be removed from the model
-0100 %   when using this approach, since checkTasks/fitTasks rely on putting specific
-0101 %   constraints for each task. The INIT algorithm will remove exchange metabolites
-0102 %   if any are present. Use importModel(file,false) to import a model with
-0103 %   exchange metabolites remaining.
-0104 %
-0105 % Usage: [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks,...
-0106 %               deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=...
-0107 %               getINITModel(refModel, tissue, celltype, hpaData, arrayData,...
-0108 %               metabolomicsData, taskFile, useScoresForTasks, printReport,...
-0109 %               taskStructure, params, paramsFT)
-0110 
-0111 if nargin<3
-0112     celltype=[];
-0113 else
-0114     celltype=char(celltype);
-0115 end
-0116 if nargin<4
-0117     hpaData=[];
-0118 end
-0119 if nargin<5
-0120     arrayData=[];
-0121 end
-0122 if nargin<6
-0123     metabolomicsData=[];
-0124 end
-0125 if nargin<7
-0126     taskFile=[];
-0127 else
-0128     taskFile=char(taskFile);
-0129 end
-0130 if nargin<8 || isempty(useScoresForTasks)
-0131     useScoresForTasks=true;
-0132 end
-0133 if nargin<9 || isempty(printReport)
-0134     printReport=true;
-0135 end
-0136 if nargin<10
-0137     taskStructure=[];
-0138 end
-0139 if nargin<11
-0140     params=[];
-0141 end
-0142 if nargin<12
-0143     paramsFT=[];
-0144 end
-0145 
-0146 %Check that the model is in the closed form
-0147 if ~isfield(refModel,'unconstrained')
-0148     EM='Exchange metabolites should normally not be removed from the model when using getINITModel. Use importModel(file,false) to import a model with exchange metabolites remaining (see the documentation for details)';
-0149     dispEM(EM);
-0150 end
-0151 
-0152 %Create the task structure if not supplied
-0153 if any(taskFile) && isempty(taskStructure)
-0154     taskStructure=parseTaskList(taskFile);
-0155 end
-0156 
-0157 
-0158 % sc-tINIT to identify confidence levels of gene expression
-0159 if ~isempty(arrayData) && isfield(arrayData,'singleCells')
-0160     if arrayData.singleCells == 1
-0161         % Check to ensure cell type is defined
-0162         if ~isfield(arrayData,'celltypes')
-0163             dispEM('arrayData must contain cell type information if sc-tINIT is to be used','false');   
-0164         end
-0165         if ~ismember(upper(celltype),upper(arrayData.celltypes))
-0166             dispEM('The cell type name does not match');   
-0167         end
-0168         
-0169         % Analyze only cell type of interest
-0170         J= strcmpi(arrayData.celltypes,celltype);
-0171         
-0172         % Analyze only genes included in the reference model
-0173         I=ismember(arrayData.genes,refModel.genes);
-0174         
-0175         % Convert expression data to population fractions
-0176         binary_levels = arrayData.levels(I,J)~=0; % Classify each gene as detected (1) or not (0) in each cell
-0177         cell_count_levels = sum(binary_levels,2); % Number of cells expressing each transcript
-0178         cell_frac_levels = cell_count_levels/size(binary_levels,2); % Number of cells expressing each transcript
-0179 
-0180         % Bin cell_frac_counts manually
-0181         x = 0:.01:1;
-0182         for(i = 1:length(x))
-0183             cell_frac_count(i) = sum(round(cell_frac_levels,2)==x(i));
-0184         end
-0185         
-0186         % Fit four beta distributions
-0187         cell_frac_count(cell_frac_count==0) = NaN; % Remove zeros from optimization
-0188         cell_frac_count(1) = NaN; % Remove non-expressed genes from optimization
-0189         x_lim = 1; % Somewhat arbitrary parameter to fit left tail of distr.
-0190         myfun = @(par) nansum((cell_frac_count(1:find(x>=x_lim,1)) - ...
-0191             abs(par(1))*betapdf(x(1:find(x>=x_lim,1)),abs(par(2)),abs(par(3))) - ...
-0192             abs(par(4))*betapdf(x(1:find(x>=x_lim,1)),abs(par(5)),abs(par(6))) - ...
-0193             abs(par(7))*betapdf(x(1:find(x>=x_lim,1)),abs(par(8)),abs(par(9))) - ...
-0194             abs(par(10))*betapdf(x(1:find(x>=x_lim,1)),abs(par(11)),abs(par(12)))).^2);
-0195         
-0196         par0 = [4,2,100,7,2,30,7,5,20,5,15,20];
-0197         opts = optimset('Display','off');
-0198         [par,f_val] = fminsearch(myfun,par0,opts);
-0199         par = abs(par);
-0200         
-0201         % Plot results
-0202         if (isfield(arrayData,'plotResults'))
-0203             if arrayData.plotResults == true
-0204                 figure(); hold on; plot(x,cell_frac_count,'ko','MarkerSize',5);
-0205                 plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))),'b-','LineWidth',1)
-0206                 plot(x,abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))),'b-','LineWidth',1)
-0207                 plot(x,abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))),'b-','LineWidth',1)
-0208                 plot(x,abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'b-','LineWidth',1)
-0209                 plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))) + ...
-0210                     abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))) + ...
-0211                     abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))) + ...
-0212                     abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'-','Color',[.5 .5 .5],'LineWidth',2)
-0213                 xlabel('Expression Probability');ylabel('# of genes');set(gca,'FontSize',14,'LineWidth',1.25);
-0214                 title('Expression prediction','FontSize',18,'FontWeight','bold')
-0215             end
-0216         end
-0217         
-0218         % Score genes based on population expression (p = .05)
-0219         exprs_cutoff_1 = find(cdf('beta',x,par(2),par(3)) >.95,1)-1; % Find index of no confidence genes
-0220         exprs_cutoff_2 = find(cdf('beta',x,par(5),par(6)) >.95,1)-1; % Find index of low confidence genes
-0221         exprs_cutoff_3 = find(cdf('beta',x,par(8),par(9)) >.95,1)-1; % Find index of low confidence genes
-0222         exprs_cutoffs = sort([exprs_cutoff_1,exprs_cutoff_2,exprs_cutoff_3]);
-0223         gene_scores = cell_frac_levels*0;
-0224         gene_scores(cell_frac_levels <= x(exprs_cutoffs(1))) = 4; % Not detected
-0225         gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(1))).*(cell_frac_levels < x(exprs_cutoffs(2))))) = 3; % Low detection
-0226         gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(2))).*(cell_frac_levels < x(exprs_cutoffs(3))))) = 2; % Medium detection
-0227         gene_scores(cell_frac_levels > x(exprs_cutoffs(3))) = 1; % High detection
-0228 
-0229         % Replace hpaData with singleCellData
-0230         if printReport==true
-0231             dispEM('Single cell data is not currently compatible with HPA data. \n         Replacing hpaData with single cell-based scoring.',false);
-0232         end
-0233         hpaData.genes = arrayData.genes;
-0234         hpaData.tissues = arrayData.tissues;
-0235         hpaData.celltypes = arrayData.celltypes;
-0236         hpaData.levels = [{'High'},{'Medium'},{'Low'},{'None'}];
-0237         hpaData.gene2Level = zeros(length(arrayData.genes),length(arrayData.celltypes));
-0238         for i = 1:length(find(J))
-0239             find_var = find(J,i);
-0240             hpaData.gene2Level(I,find_var(end)) = gene_scores;
-0241         end
-0242         
-0243         % Remove arrayData from the analysis (Might be a bad idea)
-0244         clear arrayData
-0245         arrayData=[];
-0246     end
-0247 end
-0248 
-0249 
-0250 if printReport==true
-0251     if any(celltype)
-0252         fprintf(['***Generating model for: ' tissue ' - ' celltype '\n']);
-0253     else
-0254         fprintf(['***Generating model for: ' tissue '\n']);
-0255     end
-0256     if ~isempty(hpaData)
-0257         fprintf('-Using HPA data\n');
-0258     end
-0259     if ~isempty(arrayData)
-0260         fprintf('-Using array data\n');
-0261     end
-0262     if ~isempty(metabolomicsData)
-0263         fprintf('-Using metabolomics data\n');
-0264     end
-0265     if ~isempty(taskFile) || ~isempty(taskStructure)
-0266         fprintf('-Using metabolic tasks\n');
-0267     end
-0268     fprintf('\n');
-0269     
-0270     printScores(refModel,'Reference model statistics',hpaData,arrayData,tissue,celltype);
-0271 end
-0272 
-0273 %Remove dead-end reactions to speed up the optimization and to
-0274 %differentiate between reactions removed by INIT and those that are
-0275 %dead-end
-0276 [~, deletedDeadEndRxns]=simplifyModel(refModel,true,false,true,true,true);
-0277 cModel=removeReactions(refModel,deletedDeadEndRxns,false,true);
-0278 
-0279 %Store the connected model like this to keep track of stuff
-0280 if printReport==true
-0281     printScores(cModel,'Pruned model statistics',hpaData,arrayData,tissue,celltype);
-0282 end
-0283 
-0284 %If tasks have been defined, then go through them and get essential
-0285 %reactions
-0286 if ~isempty(taskStructure)
-0287     [taskReport, essentialRxnMat]=checkTasks(cModel,[],printReport,true,true,taskStructure);
-0288     
-0289     essentialRxnsForTasks=cModel.rxns(any(essentialRxnMat,2));
-0290     
-0291     %Remove tasks that cannot be performed
-0292     taskStructure(taskReport.ok==false)=[];
-0293     if printReport==true
-0294         printScores(removeReactions(cModel,setdiff(cModel.rxns,essentialRxnsForTasks),true,true),'Reactions essential for tasks',hpaData,arrayData,tissue,celltype);
-0295     end
-0296 else
-0297     essentialRxnsForTasks={};
-0298 end
-0299 
-0300 %Score the connected model
-0301 [rxnScores, geneScores]=scoreModel(cModel,hpaData,arrayData,tissue,celltype);
-0302 
-0303 %Run the INIT algorithm. The exchange reactions that are used in the final
-0304 %reactions will be open, which doesn't fit with the last step. Therefore
-0305 %delete reactions from the original model instead of taking the output. The
-0306 %default implementation does not constrain reversible reactions to only
-0307 %carry flux in one direction. Runs without the constraints on reversibility
-0308 %and with all output allowed. This is to reduce the complexity of the
-0309 %problem.
-0310 [~, deletedRxnsInINIT, metProduction]=runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
-0311 initModel=removeReactions(cModel,deletedRxnsInINIT,true,true);
-0312 if printReport==true
-0313     printScores(initModel,'INIT model statistics',hpaData,arrayData,tissue,celltype);
-0314     printScores(removeReactions(cModel,setdiff(cModel.rxns,deletedRxnsInINIT),true,true),'Reactions deleted by INIT',hpaData,arrayData,tissue,celltype);
-0315 end
-0316 
-0317 %The full model has exchange reactions in it. fitTasks calls on fillGaps,
-0318 %which automatically removes exchange metabolites (because it assumes that
-0319 %the reactions are constrained when appropriate). In this case the
-0320 %uptakes/outputs are retrieved from the task sheet instead. To prevent
-0321 %exchange reactions being used to fill gaps, they are delete from the
-0322 %reference model here.
-0323 initModel.id='INITModel';
-0324 
-0325 %If gaps in the model should be filled using a task list
-0326 if ~isempty(taskStructure)
-0327     %Remove exchange reactions and reactions already included in the INIT
-0328     %model
-0329     refModelNoExc=removeReactions(refModel,union(initModel.rxns,getExchangeRxns(refModel)),true,true);
-0330     
-0331     %At this stage the model is fully connected and most of the genes with
-0332     %good scores should have been included. The final gap-filling should
-0333     %take the scores of the genes into account, so that "rather bad"
-0334     %reactions are preferred to "very bad" reactions. However, reactions
-0335     %with positive scores will be included even if they are not connected
-0336     %in the current formulation. Therefore, such reactions will have to be
-0337     %assigned a small negative score instead.
-0338     if useScoresForTasks==true
-0339         refRxnScores=scoreModel(refModelNoExc,hpaData,arrayData,tissue,celltype);
-0340         [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,min(refRxnScores,-0.1),taskStructure,paramsFT);
-0341     else
-0342         [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,[],taskStructure,paramsFT);
-0343     end
-0344     if printReport==true
-0345         printScores(outModel,'Functional model statistics',hpaData,arrayData,tissue,celltype);
-0346         printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,arrayData,tissue,celltype);
-0347     end
-0348     
-0349     addedRxnsForTasks=refModelNoExc.rxns(any(addedRxnMat,2));
-0350 else
-0351     outModel=initModel;
-0352     addedRxnMat=[];
-0353     addedRxnsForTasks={};
-0354 end
-0355 
-0356 %The model can now perform all the tasks defined in the task list. The
-0357 %algorithm cannot deal with gene-complexes at the moment. It is therefore
-0358 %ok to remove bad genes from a reaction (as long as at least one gene is
-0359 %kept)
-0360 model=outModel;
-0361 
-0362 [~, I]=ismember(model.genes,cModel.genes); %All should be found
-0363 %This is a little weird way to make sure that only one bad gene is included
-0364 %if there are no good ones (since all -Inf==max(-Inf))
-0365 geneScores(isinf(geneScores))=-1000+rand(sum(isinf(geneScores)),1);
-0366 
-0367 model.grRules(:)={''};
-0368 for i=1:numel(model.rxns)
-0369     ids=find(model.rxnGeneMat(i,:));
-0370     if numel(ids)>1
-0371         scores=geneScores(I(ids));
-0372         %Only keep the positive ones if possible
-0373         model.rxnGeneMat(i,ids(~(scores>0 | scores==max(scores))))=0;
-0374     end
-0375     %Rewrite the grRules to be only OR
-0376     if isfield(model,'grRules')
-0377         J=find(model.rxnGeneMat(i,:));
-0378         for j=1:numel(J)
-0379             model.grRules{i}=[model.grRules{i} '(' model.genes{J(j)} ')'];
-0380             if j<numel(J)
-0381                 model.grRules{i}=[model.grRules{i} ' or '];
-0382             end
-0383         end
-0384     end
-0385 end
-0386 
-0387 %Find all genes that are not used and delete them
-0388 I=sum(model.rxnGeneMat)==0;
-0389 model.genes(I)=[];
-0390 model.rxnGeneMat(:,I)=[];
-0391 if isfield(model,'geneShortNames')
-0392     model.geneShortNames(I)=[];
-0393 end
-0394 if isfield(model,'proteins')
-0395     model.proteins(I)=[];
-0396 end
-0397 if isfield(model,'geneMiriams')
-0398     model.geneMiriams(I)=[];
-0399 end
-0400 if isfield(model,'geneFrom')
-0401     model.geneFrom(I)=[];
-0402 end
-0403 if isfield(model,'geneComps')
-0404     model.geneComps(I)=[];
-0405 end
-0406 
-0407 %At this stage the model will contain some exchange reactions but probably
-0408 %not all (and maybe zero). This can be inconvenient, so all exchange
-0409 %reactions from the reference model are added, except for those which
-0410 %involve metabolites that are not in the model.
+0056 %
+0057 % Output:
+0058 %   model                   the resulting model structure
+0059 %   metProduction           array that indicates which of the
+0060 %                           metabolites in metabolomicsData that could be
+0061 %                           produced. Note that this is before the
+0062 %                           gap-filling process to enable defined tasks. To
+0063 %                           see which metabolites that can be produced in
+0064 %                           the final model, use canProduce.
+0065 %                           -2: metabolite name not found in model
+0066 %                           -1: metabolite found, but it could not be produced
+0067 %                           1: metabolite could be produced
+0068 %   essentialRxnsForTasks   cell array of the reactions which were
+0069 %                           essential to perform the tasks
+0070 %   addedRxnsForTasks       cell array of the reactions which were added in
+0071 %                           order to perform the tasks
+0072 %   deletedDeadEndRxns      cell array of reactions deleted because they
+0073 %                           could not carry flux (INIT requires a
+0074 %                           functional input model)
+0075 %   deletedRxnsInINIT       cell array of the reactions which were deleted by
+0076 %                           the INIT algorithm
+0077 %   taskReport              structure with the results for each task
+0078 %       id                  cell array with the id of the task
+0079 %       description         cell array with the description of the task
+0080 %       ok                  boolean array with true if the task was successful
+0081 %       essential           cell array with cell arrays of essential
+0082 %                           reactions for the task
+0083 %       gapfill             cell array of cell arrays of reactions included
+0084 %                           in the gap-filling for the task
+0085 %
+0086 %   This is the main function for automatic reconstruction of models based
+0087 %   on the (t)INIT algorithm (PLoS Comput Biol. 2012;8(5):e1002518,
+0088 %   Mol Syst Biol. 2014;10:721). Not all settings are possible using this
+0089 %   function, and you may want to call the functions scoreModel, runINIT
+0090 %   and fitTasks individually instead.
+0091 %
+0092 %   NOTE: Exchange metabolites should normally not be removed from the model
+0093 %   when using this approach, since checkTasks/fitTasks rely on putting specific
+0094 %   constraints for each task. The INIT algorithm will remove exchange metabolites
+0095 %   if any are present. Use importModel(file,false) to import a model with
+0096 %   exchange metabolites remaining.
+0097 %
+0098 % Usage: [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks,...
+0099 %               deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=...
+0100 %               getINITModel(refModel, tissue, celltype, hpaData, arrayData,...
+0101 %               metabolomicsData, taskFile, useScoresForTasks, printReport,...
+0102 %               taskStructure, params, paramsFT)
+0103 
+0104 if nargin<3
+0105     celltype=[];
+0106 else
+0107     celltype=char(celltype);
+0108 end
+0109 if nargin<4
+0110     hpaData=[];
+0111 end
+0112 if nargin<5
+0113     arrayData=[];
+0114 end
+0115 if nargin<6
+0116     metabolomicsData=[];
+0117 end
+0118 if nargin<7
+0119     taskFile=[];
+0120 else
+0121     taskFile=char(taskFile);
+0122 end
+0123 if nargin<8 || isempty(useScoresForTasks)
+0124     useScoresForTasks=true;
+0125 end
+0126 if nargin<9 || isempty(printReport)
+0127     printReport=true;
+0128 end
+0129 if nargin<10
+0130     taskStructure=[];
+0131 end
+0132 if nargin<11
+0133     params=[];
+0134 end
+0135 if nargin<12
+0136     paramsFT=[];
+0137 end
+0138 
+0139 %Check that the model is in the closed form
+0140 if ~isfield(refModel,'unconstrained')
+0141     EM='Exchange metabolites should normally not be removed from the model when using getINITModel. Use importModel(file,false) to import a model with exchange metabolites remaining (see the documentation for details)';
+0142     dispEM(EM);
+0143 end
+0144 
+0145 %Create the task structure if not supplied
+0146 if any(taskFile) && isempty(taskStructure)
+0147     taskStructure=parseTaskList(taskFile);
+0148 end
+0149 
+0150 
+0151 % sc-tINIT to identify confidence levels of gene expression
+0152 if ~isempty(arrayData) && isfield(arrayData,'singleCells')
+0153     if arrayData.singleCells == 1
+0154         % Check to ensure cell type is defined
+0155         if ~isfield(arrayData,'celltypes')
+0156             dispEM('arrayData must contain cell type information if sc-tINIT is to be used','false');   
+0157         end
+0158         if ~ismember(upper(celltype),upper(arrayData.celltypes))
+0159             dispEM('The cell type name does not match');   
+0160         end
+0161         
+0162         % Analyze only cell type of interest
+0163         J= strcmpi(arrayData.celltypes,celltype);
+0164         
+0165         % Analyze only genes included in the reference model
+0166         I=ismember(arrayData.genes,refModel.genes);
+0167         
+0168         % Convert expression data to population fractions
+0169         binary_levels = arrayData.levels(I,J)~=0; % Classify each gene as detected (1) or not (0) in each cell
+0170         cell_count_levels = sum(binary_levels,2); % Number of cells expressing each transcript
+0171         cell_frac_levels = cell_count_levels/size(binary_levels,2); % Number of cells expressing each transcript
+0172 
+0173         % Bin cell_frac_counts manually
+0174         x = 0:.01:1;
+0175         for(i = 1:length(x))
+0176             cell_frac_count(i) = sum(round(cell_frac_levels,2)==x(i));
+0177         end
+0178         
+0179         % Fit four beta distributions
+0180         cell_frac_count(cell_frac_count==0) = NaN; % Remove zeros from optimization
+0181         cell_frac_count(1) = NaN; % Remove non-expressed genes from optimization
+0182         x_lim = 1; % Somewhat arbitrary parameter to fit left tail of distr.
+0183         myfun = @(par) nansum((cell_frac_count(1:find(x>=x_lim,1)) - ...
+0184             abs(par(1))*betapdf(x(1:find(x>=x_lim,1)),abs(par(2)),abs(par(3))) - ...
+0185             abs(par(4))*betapdf(x(1:find(x>=x_lim,1)),abs(par(5)),abs(par(6))) - ...
+0186             abs(par(7))*betapdf(x(1:find(x>=x_lim,1)),abs(par(8)),abs(par(9))) - ...
+0187             abs(par(10))*betapdf(x(1:find(x>=x_lim,1)),abs(par(11)),abs(par(12)))).^2);
+0188         
+0189         par0 = [4,2,100,7,2,30,7,5,20,5,15,20];
+0190         opts = optimset('Display','off');
+0191         [par,f_val] = fminsearch(myfun,par0,opts);
+0192         par = abs(par);
+0193         
+0194         % Plot results
+0195         if (isfield(arrayData,'plotResults'))
+0196             if arrayData.plotResults == true
+0197                 figure(); hold on; plot(x,cell_frac_count,'ko','MarkerSize',5);
+0198                 plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))),'b-','LineWidth',1)
+0199                 plot(x,abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))),'b-','LineWidth',1)
+0200                 plot(x,abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))),'b-','LineWidth',1)
+0201                 plot(x,abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'b-','LineWidth',1)
+0202                 plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))) + ...
+0203                     abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))) + ...
+0204                     abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))) + ...
+0205                     abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'-','Color',[.5 .5 .5],'LineWidth',2)
+0206                 xlabel('Expression Probability');ylabel('# of genes');set(gca,'FontSize',14,'LineWidth',1.25);
+0207                 title('Expression prediction','FontSize',18,'FontWeight','bold')
+0208             end
+0209         end
+0210         
+0211         % Score genes based on population expression (p = .05)
+0212         exprs_cutoff_1 = find(cdf('beta',x,par(2),par(3)) >.95,1)-1; % Find index of no confidence genes
+0213         exprs_cutoff_2 = find(cdf('beta',x,par(5),par(6)) >.95,1)-1; % Find index of low confidence genes
+0214         exprs_cutoff_3 = find(cdf('beta',x,par(8),par(9)) >.95,1)-1; % Find index of low confidence genes
+0215         exprs_cutoffs = sort([exprs_cutoff_1,exprs_cutoff_2,exprs_cutoff_3]);
+0216         gene_scores = cell_frac_levels*0;
+0217         gene_scores(cell_frac_levels <= x(exprs_cutoffs(1))) = 4; % Not detected
+0218         gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(1))).*(cell_frac_levels < x(exprs_cutoffs(2))))) = 3; % Low detection
+0219         gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(2))).*(cell_frac_levels < x(exprs_cutoffs(3))))) = 2; % Medium detection
+0220         gene_scores(cell_frac_levels > x(exprs_cutoffs(3))) = 1; % High detection
+0221 
+0222         % Replace hpaData with singleCellData
+0223         if printReport==true
+0224             dispEM('Single cell data is not currently compatible with HPA data. \n         Replacing hpaData with single cell-based scoring.',false);
+0225         end
+0226         hpaData.genes = arrayData.genes;
+0227         hpaData.tissues = arrayData.tissues;
+0228         hpaData.celltypes = arrayData.celltypes;
+0229         hpaData.levels = [{'High'},{'Medium'},{'Low'},{'None'}];
+0230         hpaData.gene2Level = zeros(length(arrayData.genes),length(arrayData.celltypes));
+0231         for i = 1:length(find(J))
+0232             find_var = find(J,i);
+0233             hpaData.gene2Level(I,find_var(end)) = gene_scores;
+0234         end
+0235         
+0236         % Remove arrayData from the analysis (Might be a bad idea)
+0237         clear arrayData
+0238         arrayData=[];
+0239     end
+0240 end
+0241 
+0242 
+0243 if printReport==true
+0244     if any(celltype)
+0245         fprintf(['***Generating model for: ' tissue ' - ' celltype '\n']);
+0246     else
+0247         fprintf(['***Generating model for: ' tissue '\n']);
+0248     end
+0249     if ~isempty(hpaData)
+0250         fprintf('-Using HPA data\n');
+0251     end
+0252     if ~isempty(arrayData)
+0253         fprintf('-Using array data\n');
+0254     end
+0255     if ~isempty(metabolomicsData)
+0256         fprintf('-Using metabolomics data\n');
+0257     end
+0258     if ~isempty(taskFile) || ~isempty(taskStructure)
+0259         fprintf('-Using metabolic tasks\n');
+0260     end
+0261     fprintf('\n');
+0262     
+0263     printScores(refModel,'Reference model statistics',hpaData,arrayData,tissue,celltype);
+0264 end
+0265 
+0266 %Remove dead-end reactions to speed up the optimization and to
+0267 %differentiate between reactions removed by INIT and those that are
+0268 %dead-end
+0269 [~, deletedDeadEndRxns]=simplifyModel(refModel,true,false,true,true,true);
+0270 cModel=removeReactions(refModel,deletedDeadEndRxns,false,true);
+0271 
+0272 %Store the connected model like this to keep track of stuff
+0273 if printReport==true
+0274     printScores(cModel,'Pruned model statistics',hpaData,arrayData,tissue,celltype);
+0275 end
+0276 
+0277 %If tasks have been defined, then go through them and get essential
+0278 %reactions
+0279 if ~isempty(taskStructure)
+0280     [taskReport, essentialRxnMat]=checkTasks(cModel,[],printReport,true,true,taskStructure);
+0281     
+0282     essentialRxnsForTasks=cModel.rxns(any(essentialRxnMat,2));
+0283     
+0284     %Remove tasks that cannot be performed
+0285     taskStructure(taskReport.ok==false)=[];
+0286     if printReport==true
+0287         printScores(removeReactions(cModel,setdiff(cModel.rxns,essentialRxnsForTasks),true,true),'Reactions essential for tasks',hpaData,arrayData,tissue,celltype);
+0288     end
+0289 else
+0290     essentialRxnsForTasks={};
+0291 end
+0292 
+0293 %Score the connected model
+0294 [rxnScores, geneScores]=scoreModel(cModel,hpaData,arrayData,tissue,celltype);
+0295 
+0296 %Run the INIT algorithm. The exchange reactions that are used in the final
+0297 %reactions will be open, which doesn't fit with the last step. Therefore
+0298 %delete reactions from the original model instead of taking the output. The
+0299 %default implementation does not constrain reversible reactions to only
+0300 %carry flux in one direction. Runs without the constraints on reversibility
+0301 %and with all output allowed. This is to reduce the complexity of the
+0302 %problem.
+0303 [~, deletedRxnsInINIT, metProduction]=runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
+0304 initModel=removeReactions(cModel,deletedRxnsInINIT,true,true);
+0305 if printReport==true
+0306     printScores(initModel,'INIT model statistics',hpaData,arrayData,tissue,celltype);
+0307     printScores(removeReactions(cModel,setdiff(cModel.rxns,deletedRxnsInINIT),true,true),'Reactions deleted by INIT',hpaData,arrayData,tissue,celltype);
+0308 end
+0309 
+0310 %The full model has exchange reactions in it. fitTasks calls on fillGaps,
+0311 %which automatically removes exchange metabolites (because it assumes that
+0312 %the reactions are constrained when appropriate). In this case the
+0313 %uptakes/outputs are retrieved from the task sheet instead. To prevent
+0314 %exchange reactions being used to fill gaps, they are delete from the
+0315 %reference model here.
+0316 initModel.id='INITModel';
+0317 
+0318 %If gaps in the model should be filled using a task list
+0319 if ~isempty(taskStructure)
+0320     %Remove exchange reactions and reactions already included in the INIT
+0321     %model
+0322     refModelNoExc=removeReactions(refModel,union(initModel.rxns,getExchangeRxns(refModel)),true,true);
+0323     
+0324     %At this stage the model is fully connected and most of the genes with
+0325     %good scores should have been included. The final gap-filling should
+0326     %take the scores of the genes into account, so that "rather bad"
+0327     %reactions are preferred to "very bad" reactions. However, reactions
+0328     %with positive scores will be included even if they are not connected
+0329     %in the current formulation. Therefore, such reactions will have to be
+0330     %assigned a small negative score instead.
+0331     if useScoresForTasks==true
+0332         refRxnScores=scoreModel(refModelNoExc,hpaData,arrayData,tissue,celltype);
+0333         [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,min(refRxnScores,-0.1),taskStructure,paramsFT);
+0334     else
+0335         [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,[],taskStructure,paramsFT);
+0336     end
+0337     if printReport==true
+0338         printScores(outModel,'Functional model statistics',hpaData,arrayData,tissue,celltype);
+0339         printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,arrayData,tissue,celltype);
+0340     end
+0341     
+0342     addedRxnsForTasks=refModelNoExc.rxns(any(addedRxnMat,2));
+0343 else
+0344     outModel=initModel;
+0345     addedRxnMat=[];
+0346     addedRxnsForTasks={};
+0347 end
+0348 
+0349 %The model can now perform all the tasks defined in the task list. The
+0350 %algorithm cannot deal with gene-complexes at the moment. It is therefore
+0351 %ok to remove bad genes from a reaction (as long as at least one gene is
+0352 %kept)
+0353 model=outModel;
+0354 
+0355 [~, I]=ismember(model.genes,cModel.genes); %All should be found
+0356 %This is a little weird way to make sure that only one bad gene is included
+0357 %if there are no good ones (since all -Inf==max(-Inf))
+0358 geneScores(isinf(geneScores))=-1000+rand(sum(isinf(geneScores)),1);
+0359 
+0360 model.grRules(:)={''};
+0361 for i=1:numel(model.rxns)
+0362     ids=find(model.rxnGeneMat(i,:));
+0363     if numel(ids)>1
+0364         scores=geneScores(I(ids));
+0365         %Only keep the positive ones if possible
+0366         model.rxnGeneMat(i,ids(~(scores>0 | scores==max(scores))))=0;
+0367     end
+0368     %Rewrite the grRules to be only OR
+0369     if isfield(model,'grRules')
+0370         J=find(model.rxnGeneMat(i,:));
+0371         for j=1:numel(J)
+0372             model.grRules{i}=[model.grRules{i} '(' model.genes{J(j)} ')'];
+0373             if j<numel(J)
+0374                 model.grRules{i}=[model.grRules{i} ' or '];
+0375             end
+0376         end
+0377     end
+0378 end
+0379 
+0380 %Find all genes that are not used and delete them
+0381 I=sum(model.rxnGeneMat)==0;
+0382 model.genes(I)=[];
+0383 model.rxnGeneMat(:,I)=[];
+0384 if isfield(model,'geneShortNames')
+0385     model.geneShortNames(I)=[];
+0386 end
+0387 if isfield(model,'proteins')
+0388     model.proteins(I)=[];
+0389 end
+0390 if isfield(model,'geneMiriams')
+0391     model.geneMiriams(I)=[];
+0392 end
+0393 if isfield(model,'geneFrom')
+0394     model.geneFrom(I)=[];
+0395 end
+0396 if isfield(model,'geneComps')
+0397     model.geneComps(I)=[];
+0398 end
+0399 
+0400 %At this stage the model will contain some exchange reactions but probably
+0401 %not all (and maybe zero). This can be inconvenient, so all exchange
+0402 %reactions from the reference model are added, except for those which
+0403 %involve metabolites that are not in the model.
+0404 
+0405 %First delete and included exchange reactions in order to prevent the order
+0406 %from changing
+0407 model=removeReactions(model,getExchangeRxns(model));
+0408 
+0409 %Create a model with only the exchange reactions in refModel
+0410 excModel=removeReactions(refModel,setdiff(refModel.rxns,getExchangeRxns(refModel)),true,true);
 0411 
-0412 %First delete and included exchange reactions in order to prevent the order
-0413 %from changing
-0414 model=removeReactions(model,getExchangeRxns(model));
+0412 %Find the metabolites there which are not exchange metabolites and which do
+0413 %not exist in the output model
+0414 I=~ismember(excModel.mets,model.mets) & excModel.unconstrained==0;
 0415 
-0416 %Create a model with only the exchange reactions in refModel
-0417 excModel=removeReactions(refModel,setdiff(refModel.rxns,getExchangeRxns(refModel)),true,true);
-0418 
-0419 %Find the metabolites there which are not exchange metabolites and which do
-0420 %not exist in the output model
-0421 I=~ismember(excModel.mets,model.mets) & excModel.unconstrained==0;
-0422 
-0423 %Then find those reactions and delete them
-0424 [~, J]=find(excModel.S(I,:));
-0425 excModel=removeReactions(excModel,J,true,true);
-0426 
-0427 %Merge with the output model
-0428 model=mergeModels({model;excModel},'metNames');
-0429 model.id='INITModel';
-0430 model.name=['Automatically generated model for ' tissue];
-0431 if any(celltype)
-0432     model.name=[model.name ' - ' celltype];
-0433 end
-0434 
-0435 if printReport==true
-0436     printScores(model,'Final model statistics',hpaData,arrayData,tissue,celltype);
-0437 end
-0438 
-0439 %Add information about essential reactions and reactions included for
-0440 %gap-filling and return a taskReport
-0441 if ~isempty(taskStructure)
-0442     I=find(taskReport.ok); %Ignore failed tasks
-0443     for i=1:numel(I)
-0444         taskReport.essential{I(i),1}=cModel.rxns(essentialRxnMat(:,I(i)));
-0445         taskReport.gapfill{I(i),1}=refModelNoExc.rxns(addedRxnMat(:,i));
-0446     end
-0447 else
-0448     taskReport=[];
-0449 end
-0450 
-0451 %Fix grRules and reconstruct rxnGeneMat
-0452 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
-0453 model.grRules = grRules;
-0454 model.rxnGeneMat = rxnGeneMat;
-0455 end
-0456 
-0457 %This is for printing a summary of a model
-0458 function [rxnS, geneS]=printScores(model,name,hpaData,arrayData,tissue,celltype)
-0459 [a, b]=scoreModel(model,hpaData,arrayData,tissue,celltype);
-0460 rxnS=mean(a);
-0461 geneS=mean(b(~isinf(b)));
-0462 fprintf([name ':\n']);
-0463 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
-0464 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
-0465 fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
-0466 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
-0467 end
    
+0416 %Then find those reactions and delete them
+0417 [~, J]=find(excModel.S(I,:));
+0418 excModel=removeReactions(excModel,J,true,true);
+0419 
+0420 %Merge with the output model
+0421 model=mergeModels({model;excModel},'metNames');
+0422 model.id='INITModel';
+0423 model.name=['Automatically generated model for ' tissue];
+0424 if any(celltype)
+0425     model.name=[model.name ' - ' celltype];
+0426 end
+0427 
+0428 if printReport==true
+0429     printScores(model,'Final model statistics',hpaData,arrayData,tissue,celltype);
+0430 end
+0431 
+0432 %Add information about essential reactions and reactions included for
+0433 %gap-filling and return a taskReport
+0434 if ~isempty(taskStructure)
+0435     I=find(taskReport.ok); %Ignore failed tasks
+0436     for i=1:numel(I)
+0437         taskReport.essential{I(i),1}=cModel.rxns(essentialRxnMat(:,I(i)));
+0438         taskReport.gapfill{I(i),1}=refModelNoExc.rxns(addedRxnMat(:,i));
+0439     end
+0440 else
+0441     taskReport=[];
+0442 end
+0443 
+0444 %Fix grRules and reconstruct rxnGeneMat
+0445 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
+0446 model.grRules = grRules;
+0447 model.rxnGeneMat = rxnGeneMat;
+0448 end
+0449 
+0450 %This is for printing a summary of a model
+0451 function [rxnS, geneS]=printScores(model,name,hpaData,arrayData,tissue,celltype)
+0452 [a, b]=scoreModel(model,hpaData,arrayData,tissue,celltype);
+0453 rxnS=mean(a);
+0454 geneS=mean(b(~isinf(b)));
+0455 fprintf([name ':\n']);
+0456 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
+0457 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
+0458 fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
+0459 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
+0460 end
    
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file
diff --git a/doc/INIT/runINIT.html b/doc/INIT/runINIT.html
index 82146b24..e4003c24 100644
--- a/doc/INIT/runINIT.html
+++ b/doc/INIT/runINIT.html
@@ -24,7 +24,7 @@ 
    PURPOSE ^runINIT

    SYNOPSIS ^

    -
    function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops,params)
    +
    function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops)

    DESCRIPTION ^

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

    SYNOPSIS ^

    -
    function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params)
    +
    function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores)

    DESCRIPTION ^

     fillGaps
@@ -52,7 +52,6 @@ 
    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
 <li><a href=dispEM	dispEM
  • getExchangeRxns	getExchangeRxns
  • getMinNrFluxes	getMinNrFluxes
  • haveFlux	haveFlux
  • mergeModels	mergeModels
  • removeReactions	removeReactions
  • simplifyModel	simplifyModel
    • 
 This function is called by:
 
+
  • fitTasks	fitTasks
  • gapReport	gapReport
    • 
 
 
 
 
 
    SOURCE CODE ^
    
-
    0001 function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params)
+
    0001 function [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores)
 0002 % fillGaps
 0003 %   Uses template model(s) to fill gaps in a model
 0004 %
@@ -130,198 +129,197 @@ 
    SOURCE CODE ^%                       The solver will try to maximize the sum of the
 0025 %                       scores for the included reactions (optional, default
 0026 %                       is -1 for all reactions)
-0027 %   params              parameter structure as used by getMILPParams (optional)
-0028 %
-0029 %   newConnected        cell array with the reactions that could be
-0030 %                       connected. This is not calulated if
-0031 %                       useModelConstraints is true
-0032 %   cannotConnect       cell array with reactions that could not be
-0033 %                       connected. This is not calculated if
-0034 %                       useModelConstraints is true
-0035 %   addedRxns           cell array with the reactions that were added from
-0036 %                       "models"
-0037 %   newModel            the model with reactions added to fill gaps
-0038 %   exitFlag            1: optimal solution found
-0039 %                      -1: no feasible solution found
-0040 %                      -2: optimization time out
-0041 %
-0042 %   This method works by merging the model to the reference model(s) and
-0043 %   checking which reactions can carry flux. All reactions that can't
-0044 %   carry flux are removed (cannotConnect).
-0045 %   If useModelConstraints is false it then solves the MILP problem of
-0046 %   minimizing the number of active reactions from the reference models
-0047 %   that are required to have flux in all the reactions in model. This
-0048 %   requires that the input model has exchange reactions present for the
-0049 %   nutrients that are needed for its metabolism. If useModelConstraints is
-0050 %   true then the problem is to include as few reactions as possible from
-0051 %   the reference models in order to satisfy the model constraints.
-0052 %   The intended use is that the user can attempt a general gap-filling using
-0053 %   useModelConstraint=false or a more targeted gap-filling by setting
-0054 %   constraints in the model structure and then use
-0055 %   useModelConstraints=true. Say that the user want to include reactions
-0056 %   so that all biomass components can be synthesized. He/she could then
-0057 %   define a biomass equation and set the lower bound to >0. Running this
-0058 %   function with useModelConstraints=true would then give the smallest set
-0059 %   of reactions that have to be included in order for the model to produce
-0060 %   biomass.
-0061 %
-0062 % Usage: [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=...
-0063 %           fillGaps(model,models,allowNetProduction,useModelConstraints,...
-0064 %           supressWarnings,rxnScores,params)
-0065 
-0066 %If the user only supplied a single template model
-0067 if ~iscell(models)
-0068     models={models};
-0069 end
-0070 
-0071 if nargin<3
-0072     allowNetProduction=false;
-0073 end
-0074 if nargin<4
-0075     useModelConstraints=false;
-0076 end
-0077 if nargin<5
-0078     supressWarnings=false;
-0079 end
-0080 if nargin<6
-0081     rxnScores=cell(numel(models),1);
-0082     for i=1:numel(models)
-0083         rxnScores{i}=ones(numel(models{i}.rxns),1)*-1;
-0084     end
-0085 end
-0086 if isempty(rxnScores)
-0087     rxnScores=cell(numel(models),1);
-0088     for i=1:numel(models)
-0089         rxnScores{i}=ones(numel(models{i}.rxns),1)*-1;
-0090     end
-0091 end
-0092 if nargin<7
-0093     params=struct();
-0094 end
-0095 
-0096 if ~iscell(rxnScores)
-0097     rxnScores={rxnScores};
-0098 end
-0099 
-0100 models=models(:);
-0101 rxnScores=rxnScores(:);
-0102 
-0103 %Check if the original model has an unconstrained field. If so, give a
-0104 %warning
-0105 if supressWarnings==false
-0106     if isfield(model,'unconstrained')
-0107         EM='This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model still has the "unconstrained" field';
-0108         dispEM(EM,false);
-0109     else
-0110         if isempty(getExchangeRxns(model,'both'))
-0111             fprintf('NOTE: This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model does not seem to contain any such reactions.\n');
-0112         end
-0113     end
-0114 end
-0115 
-0116 %Simplify the template models to remove constrained rxns. At the same time,
-0117 %check that the id of the template models isn't the same as the model. That
-0118 %would cause an error further down
-0119 for i=1:numel(models)
-0120     models{i}.rxnScores=rxnScores{i};
-0121     models{i}=simplifyModel(models{i},false,false,true);
-0122     if strcmpi(models{i}.id,model.id)
-0123         EM='The reference model(s) cannot have the same id as the model';
-0124         dispEM(EM);
-0125     end
-0126 end
-0127 
-0128 %This is a rather ugly solution to the issue that it's a bit tricky to keep
-0129 %track of which scores belong to which reactions. This requires that
-0130 %removeReactions and mergeModels are modified to check for the new field.
-0131 model.rxnScores=zeros(numel(model.rxns),1);
-0132 
-0133 %First merge all models into one big one
-0134 allModels=mergeModels([{model};models],'metNames',true);
-0135 
-0136 %Add that net production is ok
-0137 if allowNetProduction==true
-0138     %A second column in model.b means that the b field is lower and upper
-0139     %bound on the RHS.
-0140     model.b=[model.b(:,1) inf(numel(model.mets),1)];
-0141     allModels.b=[allModels.b(:,1) inf(numel(allModels.mets),1)];
-0142 end
-0143 
-0144 if useModelConstraints==true
-0145     newConnected={};
-0146     cannotConnect={};
-0147     
-0148     %Check that the input model isn't solveable without any input
-0149     sol=solveLP(model);
-0150     if ~isempty(sol.f)
-0151         addedRxns={};
-0152         newModel=model;
-0153         exitFlag=1;
-0154         return;
-0155     end
-0156     
-0157     %Then check that the merged model is solveable
-0158     sol=solveLP(allModels);
-0159     if isempty(sol.f)
-0160         EM='There are no reactions in the template model(s) that can make the model constraints satisfied';
-0161         dispEM(EM);
-0162     end
-0163     
-0164     %Remove dead ends for speed reasons. This has to be done here and
-0165     %duplicate below because there is otherwise a risk that a reaction
-0166     %which is constrained to something relevant is removed
-0167     allModels=simplifyModel(allModels,false,false,false,true,false,false,false,[],true);
-0168     allModels.c(:)=0;
-0169 else
-0170     %Remove dead ends for speed reasons
-0171     allModels=simplifyModel(allModels,false,false,false,true,false,false,false,[],true);
-0172     allModels.c(:)=0;
-0173     
-0174     %If model constraints shouldn't be used, then determine which reactions
-0175     %to force to have flux
-0176     
-0177     %Get the reactions that can carry flux in the original model
-0178     originalFlux=haveFlux(model,1);
-0179     
-0180     %For the ones that can't carry flux, see if they can do so in the
-0181     %merged model
-0182     toCheck=intersect(allModels.rxns(strcmp(allModels.rxnFrom,model.id)),model.rxns(~originalFlux));
-0183     
-0184     %Get the ones that still cannot carry flux
-0185     I=haveFlux(allModels,1,toCheck);
-0186     
-0187     %Get the reactions that can't carry flux in the original model, but can
-0188     %in the merged one
-0189     K=toCheck(I);
-0190     
-0191     %This is a temporary thing to only look at the non-reversible rxns.
-0192     %This is because all reversible rxns can have a flux in the
-0193     %irreversible model format that is used by getMinNrFluxes
-0194     [~, I]=ismember(K,model.rxns);
-0195     K(model.rev(I)~=0)=[];
-0196     
-0197     %Constrain all reactions in the original model to have a flux
-0198     allModels.lb(ismember(allModels.rxns,K))=0.1;
-0199     
-0200     %Return stuff
-0201     newConnected=K;
-0202     cannotConnect=setdiff(model.rxns(~originalFlux),newConnected);
-0203 end
-0204 
-0205 %Then minimize for the number of fluxes used. The fixed rxns doesn't need
-0206 %to participate
-0207 templateRxns=find(~strcmp(allModels.rxnFrom,model.id));
-0208 [~, J, exitFlag]=getMinNrFluxes(allModels,templateRxns,params,allModels.rxnScores(templateRxns));
-0209 
-0210 %Remove everything except for the added ones
-0211 I=true(numel(allModels.rxns),1);
-0212 I(templateRxns(J))=false;
-0213 addedModel=removeReactions(allModels,I,true,true,true);
-0214 
-0215 newModel=mergeModels({model;addedModel},'metNames',true);
-0216 addedRxns=setdiff(newModel.rxns,model.rxns);
-0217 newModel=rmfield(newModel,'rxnScores');
-0218 end
    
+0027 %
+0028 %   newConnected        cell array with the reactions that could be
+0029 %                       connected. This is not calulated if
+0030 %                       useModelConstraints is true
+0031 %   cannotConnect       cell array with reactions that could not be
+0032 %                       connected. This is not calculated if
+0033 %                       useModelConstraints is true
+0034 %   addedRxns           cell array with the reactions that were added from
+0035 %                       "models"
+0036 %   newModel            the model with reactions added to fill gaps
+0037 %   exitFlag            1: optimal solution found
+0038 %                      -1: no feasible solution found
+0039 %                      -2: optimization time out
+0040 %
+0041 %   This method works by merging the model to the reference model(s) and
+0042 %   checking which reactions can carry flux. All reactions that can't
+0043 %   carry flux are removed (cannotConnect).
+0044 %   If useModelConstraints is false it then solves the MILP problem of
+0045 %   minimizing the number of active reactions from the reference models
+0046 %   that are required to have flux in all the reactions in model. This
+0047 %   requires that the input model has exchange reactions present for the
+0048 %   nutrients that are needed for its metabolism. If useModelConstraints is
+0049 %   true then the problem is to include as few reactions as possible from
+0050 %   the reference models in order to satisfy the model constraints.
+0051 %   The intended use is that the user can attempt a general gap-filling using
+0052 %   useModelConstraint=false or a more targeted gap-filling by setting
+0053 %   constraints in the model structure and then use
+0054 %   useModelConstraints=true. Say that the user want to include reactions
+0055 %   so that all biomass components can be synthesized. He/she could then
+0056 %   define a biomass equation and set the lower bound to >0. Running this
+0057 %   function with useModelConstraints=true would then give the smallest set
+0058 %   of reactions that have to be included in order for the model to produce
+0059 %   biomass.
+0060 %
+0061 % Usage: [newConnected, cannotConnect, addedRxns, newModel, exitFlag]=...
+0062 %           fillGaps(model,models,allowNetProduction,useModelConstraints,...
+0063 %           supressWarnings,rxnScores,params)
+0064 
+0065 %If the user only supplied a single template model
+0066 if ~iscell(models)
+0067     models={models};
+0068 end
+0069 
+0070 if nargin<3
+0071     allowNetProduction=false;
+0072 end
+0073 if nargin<4
+0074     useModelConstraints=false;
+0075 end
+0076 if nargin<5
+0077     supressWarnings=false;
+0078 end
+0079 if nargin<6
+0080     rxnScores=cell(numel(models),1);
+0081     for i=1:numel(models)
+0082         rxnScores{i}=ones(numel(models{i}.rxns),1)*-1;
+0083     end
+0084 end
+0085 if isempty(rxnScores)
+0086     rxnScores=cell(numel(models),1);
+0087     for i=1:numel(models)
+0088         rxnScores{i}=ones(numel(models{i}.rxns),1)*-1;
+0089     end
+0090 end
+0091 if nargin<7
+0092     params=struct();
+0093 end
+0094 
+0095 if ~iscell(rxnScores)
+0096     rxnScores={rxnScores};
+0097 end
+0098 
+0099 models=models(:);
+0100 rxnScores=rxnScores(:);
+0101 
+0102 %Check if the original model has an unconstrained field. If so, give a
+0103 %warning
+0104 if supressWarnings==false
+0105     if isfield(model,'unconstrained')
+0106         EM='This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model still has the "unconstrained" field';
+0107         dispEM(EM,false);
+0108     else
+0109         if isempty(getExchangeRxns(model,'both'))
+0110             fprintf('NOTE: This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model does not seem to contain any such reactions.\n');
+0111         end
+0112     end
+0113 end
+0114 
+0115 %Simplify the template models to remove constrained rxns. At the same time,
+0116 %check that the id of the template models isn't the same as the model. That
+0117 %would cause an error further down
+0118 for i=1:numel(models)
+0119     models{i}.rxnScores=rxnScores{i};
+0120     models{i}=simplifyModel(models{i},false,false,true);
+0121     if strcmpi(models{i}.id,model.id)
+0122         EM='The reference model(s) cannot have the same id as the model';
+0123         dispEM(EM);
+0124     end
+0125 end
+0126 
+0127 %This is a rather ugly solution to the issue that it's a bit tricky to keep
+0128 %track of which scores belong to which reactions. This requires that
+0129 %removeReactions and mergeModels are modified to check for the new field.
+0130 model.rxnScores=zeros(numel(model.rxns),1);
+0131 
+0132 %First merge all models into one big one
+0133 allModels=mergeModels([{model};models],'metNames',true);
+0134 
+0135 %Add that net production is ok
+0136 if allowNetProduction==true
+0137     %A second column in model.b means that the b field is lower and upper
+0138     %bound on the RHS.
+0139     model.b=[model.b(:,1) inf(numel(model.mets),1)];
+0140     allModels.b=[allModels.b(:,1) inf(numel(allModels.mets),1)];
+0141 end
+0142 
+0143 if useModelConstraints==true
+0144     newConnected={};
+0145     cannotConnect={};
+0146     
+0147     %Check that the input model isn't solveable without any input
+0148     sol=solveLP(model);
+0149     if ~isempty(sol.f)
+0150         addedRxns={};
+0151         newModel=model;
+0152         exitFlag=1;
+0153         return;
+0154     end
+0155     
+0156     %Then check that the merged model is solveable
+0157     sol=solveLP(allModels);
+0158     if isempty(sol.f)
+0159         EM='There are no reactions in the template model(s) that can make the model constraints satisfied';
+0160         dispEM(EM);
+0161     end
+0162     
+0163     %Remove dead ends for speed reasons. This has to be done here and
+0164     %duplicate below because there is otherwise a risk that a reaction
+0165     %which is constrained to something relevant is removed
+0166     allModels=simplifyModel(allModels,false,false,false,true,false,false,false,[],true);
+0167     allModels.c(:)=0;
+0168 else
+0169     %Remove dead ends for speed reasons
+0170     allModels=simplifyModel(allModels,false,false,false,true,false,false,false,[],true);
+0171     allModels.c(:)=0;
+0172     
+0173     %If model constraints shouldn't be used, then determine which reactions
+0174     %to force to have flux
+0175     
+0176     %Get the reactions that can carry flux in the original model
+0177     originalFlux=haveFlux(model,1);
+0178     
+0179     %For the ones that can't carry flux, see if they can do so in the
+0180     %merged model
+0181     toCheck=intersect(allModels.rxns(strcmp(allModels.rxnFrom,model.id)),model.rxns(~originalFlux));
+0182     
+0183     %Get the ones that still cannot carry flux
+0184     I=haveFlux(allModels,1,toCheck);
+0185     
+0186     %Get the reactions that can't carry flux in the original model, but can
+0187     %in the merged one
+0188     K=toCheck(I);
+0189     
+0190     %This is a temporary thing to only look at the non-reversible rxns.
+0191     %This is because all reversible rxns can have a flux in the
+0192     %irreversible model format that is used by getMinNrFluxes
+0193     [~, I]=ismember(K,model.rxns);
+0194     K(model.rev(I)~=0)=[];
+0195     
+0196     %Constrain all reactions in the original model to have a flux
+0197     allModels.lb(ismember(allModels.rxns,K))=0.1;
+0198     
+0199     %Return stuff
+0200     newConnected=K;
+0201     cannotConnect=setdiff(model.rxns(~originalFlux),newConnected);
+0202 end
+0203 
+0204 %Then minimize for the number of fluxes used. The fixed rxns doesn't need
+0205 %to participate
+0206 templateRxns=find(~strcmp(allModels.rxnFrom,model.id));
+0207 [~, J, exitFlag]=getMinNrFluxes(allModels,templateRxns,params,allModels.rxnScores(templateRxns));
+0208 
+0209 %Remove everything except for the added ones
+0210 I=true(numel(allModels.rxns),1);
+0211 I(templateRxns(J))=false;
+0212 addedModel=removeReactions(allModels,I,true,true,true);
+0213 
+0214 newModel=mergeModels({model;addedModel},'metNames',true);
+0215 addedRxns=setdiff(newModel.rxns,model.rxns);
+0216 newModel=rmfield(newModel,'rxnScores');
+0217 end
    
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file
diff --git a/doc/core/fitTasks.html b/doc/core/fitTasks.html
index 1090f858..fa6bdbcb 100644
--- a/doc/core/fitTasks.html
+++ b/doc/core/fitTasks.html
@@ -24,7 +24,7 @@ 
    PURPOSE ^fitTasks

    SYNOPSIS ^

    -
    function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params)
    +
    function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure)

    DESCRIPTION ^

     fitTasks
@@ -45,7 +45,6 @@ 
    DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
    
 This function calls:
 
+
  • addMets	addMets
  • addRxns	addRxns
  • dispEM	dispEM
  • fillGaps	fillGaps
  • mergeModels	mergeModels
  • parseTaskList	parseTaskList
  • printFluxes	printFluxes
  • removeReactions	removeReactions
  • setParam	setParam
    • 
 This function is called by:
 
@@ -78,7 +77,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 
 
 <h2><a name=SOURCE CODE ^
    
-
    0001 function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure,params)
+
    0001 function [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,rxnScores,taskStructure)
 0002 % fitTasks
 0003 %   Fills gaps in a model by including reactions from a reference model,
 0004 %   so that the resulting model can perform all the tasks in a task list
@@ -97,318 +96,317 @@ 
    SOURCE CODE ^%                   reactions (optional, default is -1 for all reactions)
 0018 %   taskStructure   structure with the tasks, as from parseTaskList. If
 0019 %                   this is supplied then inputFile is ignored (optional)
-0020 %   params          parameter structure as used by getMILPParams (optional)
+0020 %
 0021 %
-0022 %
-0023 %   Output:
-0024 %   outModel        model structure with reactions added to perform the
-0025 %                   tasks
-0026 %   addedRxns       MxN matrix with the added reactions (M) from refModel
-0027 %                   for each task (N). An element is true if the corresponding
-0028 %                   reaction is added in the corresponding task.
-0029 %                   Failed tasks and SHOULD FAIL tasks are ignored
-0030 %
-0031 %   This function fills gaps in a model by using a reference model, so
-0032 %   that the resulting model can perform a list of metabolic tasks. The
-0033 %   gap-filling is done in a task-by-task manner, rather than solving for
-0034 %   all tasks at once. This means that the order of the tasks could influence
-0035 %   the result.
-0036 %
-0037 % Usage: [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
-0038 %           rxnScores,taskStructure,params)
-0039 
-0040 if nargin<4
-0041     printOutput=true;
-0042 end
-0043 if nargin<5
-0044     rxnScores=ones(numel(refModel.rxns),1)*-1;
-0045 end
-0046 if isempty(rxnScores)
-0047     rxnScores=ones(numel(refModel.rxns),1)*-1;
-0048 end
-0049 if nargin<6
-0050     taskStructure=[];
-0051 end
-0052 if nargin<7
-0053     params=[];
-0054 end
-0055 
-0056 if isempty(taskStructure) && ~isfile(inputFile)
-0057     error('Task file %s cannot be found',string(inputFile));
-0058 end
-0059 
-0060 if strcmpi(model.id,refModel.id)
-0061     fprintf('NOTE: The model and reference model have the same IDs. The ID for the reference model was set to "refModel" in order to keep track of the origin of reactions.\n');
-0062     refModel.id='refModel';
-0063 end
-0064 
-0065 if any(rxnScores>=0)
-0066     EM='Only negative values are allowed in rxnScores';
-0067     dispEM(EM);
-0068 end
-0069 
-0070 %Prepare the input models a little
-0071 model.b=zeros(numel(model.mets),2);
-0072 modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
-0073 %This is the mets in the reference model. Used if the tasks involve
-0074 %metabolites that doesn't exist in the model
-0075 largeModelMets=upper(strcat(refModel.metNames,'[',refModel.comps(refModel.metComps),']'));
-0076 
-0077 if ~isfield(model,'unconstrained')
-0078     EM='Exchange metabolites should normally not be removed from the model when using checkTasks. Inputs and outputs are defined in the task file instead. Use importModel(file,false) to import a model with exchange metabolites remaining';
-0079     dispEM(EM,false);
-0080 end
-0081 
-0082 if isempty(taskStructure)
-0083     taskStructure=parseTaskList(inputFile);
-0084 end
-0085 
-0086 tModel=model;
-0087 addedRxns=false(numel(refModel.rxns),numel(taskStructure));
-0088 supressWarnings=false;
-0089 nAdded=0;
-0090 for i=1:numel(taskStructure)
-0091     if ~taskStructure(i).shouldFail
-0092         %Set the inputs
-0093         if ~isempty(taskStructure(i).inputs)
-0094             [I, J]=ismember(upper(taskStructure(i).inputs),modelMets);
-0095             K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
-0096             L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
-0097             %Check that all metabolites are either real metabolites or
-0098             %ALLMETS/ALLMETSIN
-0099             goodMets=I|K|L;
-0100             if ~all(goodMets)
-0101                 %Not all of the inputs could be found in the small model.
-0102                 %Check if they exist in the large model
-0103                 [found, metMatch]=ismember(upper(taskStructure(i).inputs(~goodMets)),largeModelMets);
-0104                 if ~all(found)
-0105                     EM=['Could not find all inputs in "[' taskStructure(i).id '] ' taskStructure(i).description '" in either model'];
-0106                     disp(EM);
-0107                 else
-0108                     %Otherwise add them to the model
-0109                     met.metNames=refModel.metNames(metMatch);
-0110                     met.compartments=refModel.comps(refModel.metComps(metMatch));
-0111                     
-0112                     %Add the metabolite both to the base model and the
-0113                     %model used in the current task
-0114                     model=addMets(model,met);
-0115                     tModel=addMets(tModel,met);
-0116                     modelMets=[modelMets;upper(taskStructure(i).inputs(~goodMets))];
-0117                 end
-0118                 
-0119                 %By now the indexes might be getting a bit confusing, but
-0120                 %this is to update the indexes of the "real" metabolites to
-0121                 %point to the newly added ones
-0122                 I(~goodMets)=true; %All the bad ones are fixed at this stage
-0123                 J(~goodMets)=numel(modelMets)-numel(metMatch)+1:numel(modelMets);
-0124             end
-0125             if numel(J(I))~=numel(unique(J(I)))
-0126                 EM=['The constraints on some input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
-0127                 dispEM(EM);
-0128             end
-0129             %If all metabolites should be added
-0130             if any(K)
-0131                 %Check if ALLMETS is the first metabolite. Otherwise print
-0132                 %a warning since it will write over any other constraints
-0133                 %that are set
-0134                 if K(1)==0
-0135                     EM=['ALLMETS is used as an input in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
-0136                     dispEM(EM,false);
-0137                 end
-0138                 %Use the first match of ALLMETS. There should only be one,
-0139                 %but still..
-0140                 tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
-0141             end
-0142             %If metabolites in a specific compartment should be used
-0143             if any(L)
-0144                 L=find(L);
-0145                 for j=1:numel(L)
-0146                     %The compartment defined
-0147                     compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
-0148                     %Check if it exists in the model
-0149                     C=find(ismember(upper(model.comps),compartment));
-0150                     if any(C)
-0151                         %Match to metabolites
-0152                         tModel.b(model.metComps==C,1)=taskStructure(i).UBin(L(j))*-1;
-0153                     else
-0154                         EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
-0155                         dispEM(EM);
-0156                     end
-0157                 end
-0158             end
-0159             %Then add the normal constraints
-0160             if any(J(I))
-0161                 tModel.b(J(I),1)=taskStructure(i).UBin(I)*-1;
-0162                 tModel.b(J(I),2)=taskStructure(i).LBin(I)*-1;
-0163             end
-0164         end
-0165         %Set the outputs
-0166         if ~isempty(taskStructure(i).outputs)
-0167             [I, J]=ismember(upper(taskStructure(i).outputs),modelMets);
-0168             K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
-0169             L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
-0170             %Check that all metabolites are either real metabolites or
-0171             %ALLMETS/ALLMETSIN
-0172             goodMets=I|K|L;
-0173             if ~all(goodMets)
-0174                 %Not all of the outputs could be found in the small model.
-0175                 %Check if they exist in the large model
-0176                 [found, metMatch]=ismember(upper(taskStructure(i).outputs(~goodMets)),largeModelMets);
-0177                 if ~all(found)
-0178                     EM=['Could not find all outputs in "[' taskStructure(i).id '] ' taskStructure(i).description '" in either model'];
-0179                     dispEM(EM);
-0180                 else
-0181                     %Otherwise add them to the model
-0182                     met.metNames=refModel.metNames(metMatch);
-0183                     met.compartments=refModel.comps(refModel.metComps(metMatch));
-0184                     
-0185                     %Add the metabolite both to the base model and the
-0186                     %model used in the current task
-0187                     model=addMets(model,met);
-0188                     tModel=addMets(tModel,met);
-0189                     modelMets=[modelMets;upper(taskStructure(i).outputs(~goodMets))];
-0190                 end
-0191                 
-0192                 %By now the indexes might be getting a bit confusing, but
-0193                 %this is to update the indexes of the "real" metabolites to
-0194                 %point to the newly added ones
-0195                 I(~goodMets)=true; %All the bad ones are fixed at this stage
-0196                 J(~goodMets)=numel(modelMets)-numel(metMatch)+1:numel(modelMets);
-0197             end
-0198             if numel(J(I))~=numel(unique(J(I)))
-0199                 EM=['The constraints on some output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
-0200                 dispEM(EM);
-0201             end
-0202             %If all metabolites should be added
-0203             if any(K)
-0204                 %Check if ALLMETS is the first metabolite. Otherwise print
-0205                 %a warning since it will write over any other constraints
-0206                 %that are set
-0207                 if K(1)==0
-0208                     EM=['ALLMETS is used as an output in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
-0209                     dispEM(EM,false);
-0210                 end
-0211                 %Use the first match of ALLMETS. There should only be one,
-0212                 %but still..
-0213                 tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
-0214             end
-0215             %If metabolites in a specific compartment should be used
-0216             if any(L)
-0217                 L=find(L);
-0218                 for j=1:numel(L)
-0219                     %The compartment defined
-0220                     compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
-0221                     %Check if it exists in the model
-0222                     C=find(ismember(upper(model.comps),compartment));
-0223                     if any(C)
-0224                         %Match to metabolites
-0225                         tModel.b(model.metComps==C,2)=taskStructure(i).UBout(L(j));
-0226                     else
-0227                         EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
-0228                         dispEM(EM);
-0229                     end
-0230                 end
-0231             end
-0232             %Then add the normal constraints
-0233             if any(J(I))
-0234                 %Verify that IN and OUT bounds are consistent. Cannot require
-0235                 %that a metabolite is simultaneously input AND output at some
-0236                 %nonzero flux.
-0237                 J = J(I);
-0238                 I = find(I);  % otherwise indexing becomes confusing
-0239                 nonzero_LBin = tModel.b(J,2) < 0;
-0240                 nonzero_LBout = taskStructure(i).LBout(I) > 0;
-0241                 if any(nonzero_LBin & nonzero_LBout)
-0242                     EM=['The IN LB and OUT LB in "[' taskStructure(i).id '] ' taskStructure(i).description '" cannot be nonzero for the same metabolite'];
-0243                     dispEM(EM);
-0244                 end
-0245                 tModel.b(J(nonzero_LBout),1)=taskStructure(i).LBout(I(nonzero_LBout));
-0246                 tModel.b(J,2)=taskStructure(i).UBout(I);
-0247             end
-0248         end
-0249         
-0250         %Add new rxns
-0251         if ~isempty(taskStructure(i).equations)
-0252             rxn.equations=taskStructure(i).equations;
-0253             rxn.lb=taskStructure(i).LBequ;
-0254             rxn.ub=taskStructure(i).UBequ;
-0255             rxn.rxns=strcat({'TEMPORARY_'},num2str((1:numel(taskStructure(i).equations))'));
-0256             tModel=addRxns(tModel,rxn,3);
-0257         end
-0258         %Add changed bounds
-0259         if ~isempty(taskStructure(i).changed)
-0260             tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
-0261             tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
-0262         end
-0263         
-0264         %Solve and print. Display a warning if the problem is not solveable
-0265         sol=solveLP(tModel);
-0266         if isempty(sol.x)
-0267             %Only do gap-filling if it cannot be solved
-0268             failed=false;
-0269             try
-0270                 [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params);
-0271                 if exitFlag==-2
-0272                     EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n'];
-0273                     dispEM(EM,false);
-0274                 end
-0275             catch
-0276                 EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" could not be performed for any set of reactions\n'];
-0277                 dispEM(EM,false);
-0278                 failed=true;
-0279             end
-0280             if failed==false
-0281                 if ~isempty(newRxns)
-0282                     nAdded=nAdded+numel(newRxns);
-0283                     
-0284                     %Add the reactions to the base model. It is not correct
-0285                     %to use newModel directly, as it may contain
-0286                     %reactions/constraints that are specific to this task
-0287                     model=mergeModels({model,removeReactions(newModel,setdiff(newModel.rxns,newRxns),true,true)},'metNames',true);
-0288                     
-0289                     %Keep track of the added reactions
-0290                     addedRxns(ismember(refModel.rxns,newRxns),i)=true;
-0291                 end
-0292                 if printOutput==true
-0293                     fprintf(['[' taskStructure(i).id '] ' taskStructure(i).description ': Added ' num2str(numel(newRxns)) ' reaction(s), ' num2str(nAdded) ' reactions added in total\n']);
-0294                 end
-0295             end
-0296         else
-0297             if printOutput==true
-0298                 fprintf(['[' taskStructure(i).id '] ' taskStructure(i).description ': Added 0 reaction(s), ' num2str(nAdded) ' reactions added in total\n']);
-0299             end
-0300         end
-0301         supressWarnings=true;
-0302         
-0303         %Print the output if chosen
-0304         if taskStructure(i).printFluxes && printOutput
-0305             if ~isempty(sol.x)
-0306                 sol=solveLP(tModel,1);
-0307                 printFluxes(tModel,sol.x,false,10^-5,[],'%rxnID (%eqn):%flux\n');
-0308                 fprintf('\n');
-0309             else
-0310                 %If the problem wasn't solveable then the gap-filled model
-0311                 %should be used
-0312                 if failed==false
-0313                     sol=solveLP(newModel,1);
-0314                     printFluxes(newModel,sol.x,false,10^-5,[],'%rxnID (%eqn):%flux\n');
-0315                     fprintf('\n');
-0316                 end
-0317             end
-0318         end
-0319         
-0320         tModel=model;
-0321         %Since new mets are added by merging the new reactions and not only
-0322         %from the task sheet
-0323         modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
-0324     else
-0325         EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" is set as SHOULD FAIL. Such tasks cannot be modelled using this approach and the task is therefore ignored\n'];
-0326         dispEM(EM,false);
-0327     end
-0328 end
-0329 model.b(:,2) = [];  % resume field b
-0330 outModel=model;
-0331 end
    
+0022 %   Output:
+0023 %   outModel        model structure with reactions added to perform the
+0024 %                   tasks
+0025 %   addedRxns       MxN matrix with the added reactions (M) from refModel
+0026 %                   for each task (N). An element is true if the corresponding
+0027 %                   reaction is added in the corresponding task.
+0028 %                   Failed tasks and SHOULD FAIL tasks are ignored
+0029 %
+0030 %   This function fills gaps in a model by using a reference model, so
+0031 %   that the resulting model can perform a list of metabolic tasks. The
+0032 %   gap-filling is done in a task-by-task manner, rather than solving for
+0033 %   all tasks at once. This means that the order of the tasks could influence
+0034 %   the result.
+0035 %
+0036 % Usage: [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
+0037 %           rxnScores,taskStructure,params)
+0038 
+0039 if nargin<4
+0040     printOutput=true;
+0041 end
+0042 if nargin<5
+0043     rxnScores=ones(numel(refModel.rxns),1)*-1;
+0044 end
+0045 if isempty(rxnScores)
+0046     rxnScores=ones(numel(refModel.rxns),1)*-1;
+0047 end
+0048 if nargin<6
+0049     taskStructure=[];
+0050 end
+0051 if nargin<7
+0052     params=[];
+0053 end
+0054 
+0055 if isempty(taskStructure) && ~isfile(inputFile)
+0056     error('Task file %s cannot be found',string(inputFile));
+0057 end
+0058 
+0059 if strcmpi(model.id,refModel.id)
+0060     fprintf('NOTE: The model and reference model have the same IDs. The ID for the reference model was set to "refModel" in order to keep track of the origin of reactions.\n');
+0061     refModel.id='refModel';
+0062 end
+0063 
+0064 if any(rxnScores>=0)
+0065     EM='Only negative values are allowed in rxnScores';
+0066     dispEM(EM);
+0067 end
+0068 
+0069 %Prepare the input models a little
+0070 model.b=zeros(numel(model.mets),2);
+0071 modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
+0072 %This is the mets in the reference model. Used if the tasks involve
+0073 %metabolites that doesn't exist in the model
+0074 largeModelMets=upper(strcat(refModel.metNames,'[',refModel.comps(refModel.metComps),']'));
+0075 
+0076 if ~isfield(model,'unconstrained')
+0077     EM='Exchange metabolites should normally not be removed from the model when using checkTasks. Inputs and outputs are defined in the task file instead. Use importModel(file,false) to import a model with exchange metabolites remaining';
+0078     dispEM(EM,false);
+0079 end
+0080 
+0081 if isempty(taskStructure)
+0082     taskStructure=parseTaskList(inputFile);
+0083 end
+0084 
+0085 tModel=model;
+0086 addedRxns=false(numel(refModel.rxns),numel(taskStructure));
+0087 supressWarnings=false;
+0088 nAdded=0;
+0089 for i=1:numel(taskStructure)
+0090     if ~taskStructure(i).shouldFail
+0091         %Set the inputs
+0092         if ~isempty(taskStructure(i).inputs)
+0093             [I, J]=ismember(upper(taskStructure(i).inputs),modelMets);
+0094             K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
+0095             L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
+0096             %Check that all metabolites are either real metabolites or
+0097             %ALLMETS/ALLMETSIN
+0098             goodMets=I|K|L;
+0099             if ~all(goodMets)
+0100                 %Not all of the inputs could be found in the small model.
+0101                 %Check if they exist in the large model
+0102                 [found, metMatch]=ismember(upper(taskStructure(i).inputs(~goodMets)),largeModelMets);
+0103                 if ~all(found)
+0104                     EM=['Could not find all inputs in "[' taskStructure(i).id '] ' taskStructure(i).description '" in either model'];
+0105                     disp(EM);
+0106                 else
+0107                     %Otherwise add them to the model
+0108                     met.metNames=refModel.metNames(metMatch);
+0109                     met.compartments=refModel.comps(refModel.metComps(metMatch));
+0110                     
+0111                     %Add the metabolite both to the base model and the
+0112                     %model used in the current task
+0113                     model=addMets(model,met);
+0114                     tModel=addMets(tModel,met);
+0115                     modelMets=[modelMets;upper(taskStructure(i).inputs(~goodMets))];
+0116                 end
+0117                 
+0118                 %By now the indexes might be getting a bit confusing, but
+0119                 %this is to update the indexes of the "real" metabolites to
+0120                 %point to the newly added ones
+0121                 I(~goodMets)=true; %All the bad ones are fixed at this stage
+0122                 J(~goodMets)=numel(modelMets)-numel(metMatch)+1:numel(modelMets);
+0123             end
+0124             if numel(J(I))~=numel(unique(J(I)))
+0125                 EM=['The constraints on some input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
+0126                 dispEM(EM);
+0127             end
+0128             %If all metabolites should be added
+0129             if any(K)
+0130                 %Check if ALLMETS is the first metabolite. Otherwise print
+0131                 %a warning since it will write over any other constraints
+0132                 %that are set
+0133                 if K(1)==0
+0134                     EM=['ALLMETS is used as an input in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
+0135                     dispEM(EM,false);
+0136                 end
+0137                 %Use the first match of ALLMETS. There should only be one,
+0138                 %but still..
+0139                 tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
+0140             end
+0141             %If metabolites in a specific compartment should be used
+0142             if any(L)
+0143                 L=find(L);
+0144                 for j=1:numel(L)
+0145                     %The compartment defined
+0146                     compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
+0147                     %Check if it exists in the model
+0148                     C=find(ismember(upper(model.comps),compartment));
+0149                     if any(C)
+0150                         %Match to metabolites
+0151                         tModel.b(model.metComps==C,1)=taskStructure(i).UBin(L(j))*-1;
+0152                     else
+0153                         EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
+0154                         dispEM(EM);
+0155                     end
+0156                 end
+0157             end
+0158             %Then add the normal constraints
+0159             if any(J(I))
+0160                 tModel.b(J(I),1)=taskStructure(i).UBin(I)*-1;
+0161                 tModel.b(J(I),2)=taskStructure(i).LBin(I)*-1;
+0162             end
+0163         end
+0164         %Set the outputs
+0165         if ~isempty(taskStructure(i).outputs)
+0166             [I, J]=ismember(upper(taskStructure(i).outputs),modelMets);
+0167             K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
+0168             L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
+0169             %Check that all metabolites are either real metabolites or
+0170             %ALLMETS/ALLMETSIN
+0171             goodMets=I|K|L;
+0172             if ~all(goodMets)
+0173                 %Not all of the outputs could be found in the small model.
+0174                 %Check if they exist in the large model
+0175                 [found, metMatch]=ismember(upper(taskStructure(i).outputs(~goodMets)),largeModelMets);
+0176                 if ~all(found)
+0177                     EM=['Could not find all outputs in "[' taskStructure(i).id '] ' taskStructure(i).description '" in either model'];
+0178                     dispEM(EM);
+0179                 else
+0180                     %Otherwise add them to the model
+0181                     met.metNames=refModel.metNames(metMatch);
+0182                     met.compartments=refModel.comps(refModel.metComps(metMatch));
+0183                     
+0184                     %Add the metabolite both to the base model and the
+0185                     %model used in the current task
+0186                     model=addMets(model,met);
+0187                     tModel=addMets(tModel,met);
+0188                     modelMets=[modelMets;upper(taskStructure(i).outputs(~goodMets))];
+0189                 end
+0190                 
+0191                 %By now the indexes might be getting a bit confusing, but
+0192                 %this is to update the indexes of the "real" metabolites to
+0193                 %point to the newly added ones
+0194                 I(~goodMets)=true; %All the bad ones are fixed at this stage
+0195                 J(~goodMets)=numel(modelMets)-numel(metMatch)+1:numel(modelMets);
+0196             end
+0197             if numel(J(I))~=numel(unique(J(I)))
+0198                 EM=['The constraints on some output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
+0199                 dispEM(EM);
+0200             end
+0201             %If all metabolites should be added
+0202             if any(K)
+0203                 %Check if ALLMETS is the first metabolite. Otherwise print
+0204                 %a warning since it will write over any other constraints
+0205                 %that are set
+0206                 if K(1)==0
+0207                     EM=['ALLMETS is used as an output in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
+0208                     dispEM(EM,false);
+0209                 end
+0210                 %Use the first match of ALLMETS. There should only be one,
+0211                 %but still..
+0212                 tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
+0213             end
+0214             %If metabolites in a specific compartment should be used
+0215             if any(L)
+0216                 L=find(L);
+0217                 for j=1:numel(L)
+0218                     %The compartment defined
+0219                     compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
+0220                     %Check if it exists in the model
+0221                     C=find(ismember(upper(model.comps),compartment));
+0222                     if any(C)
+0223                         %Match to metabolites
+0224                         tModel.b(model.metComps==C,2)=taskStructure(i).UBout(L(j));
+0225                     else
+0226                         EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
+0227                         dispEM(EM);
+0228                     end
+0229                 end
+0230             end
+0231             %Then add the normal constraints
+0232             if any(J(I))
+0233                 %Verify that IN and OUT bounds are consistent. Cannot require
+0234                 %that a metabolite is simultaneously input AND output at some
+0235                 %nonzero flux.
+0236                 J = J(I);
+0237                 I = find(I);  % otherwise indexing becomes confusing
+0238                 nonzero_LBin = tModel.b(J,2) < 0;
+0239                 nonzero_LBout = taskStructure(i).LBout(I) > 0;
+0240                 if any(nonzero_LBin & nonzero_LBout)
+0241                     EM=['The IN LB and OUT LB in "[' taskStructure(i).id '] ' taskStructure(i).description '" cannot be nonzero for the same metabolite'];
+0242                     dispEM(EM);
+0243                 end
+0244                 tModel.b(J(nonzero_LBout),1)=taskStructure(i).LBout(I(nonzero_LBout));
+0245                 tModel.b(J,2)=taskStructure(i).UBout(I);
+0246             end
+0247         end
+0248         
+0249         %Add new rxns
+0250         if ~isempty(taskStructure(i).equations)
+0251             rxn.equations=taskStructure(i).equations;
+0252             rxn.lb=taskStructure(i).LBequ;
+0253             rxn.ub=taskStructure(i).UBequ;
+0254             rxn.rxns=strcat({'TEMPORARY_'},num2str((1:numel(taskStructure(i).equations))'));
+0255             tModel=addRxns(tModel,rxn,3);
+0256         end
+0257         %Add changed bounds
+0258         if ~isempty(taskStructure(i).changed)
+0259             tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
+0260             tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
+0261         end
+0262         
+0263         %Solve and print. Display a warning if the problem is not solveable
+0264         sol=solveLP(tModel);
+0265         if isempty(sol.x)
+0266             %Only do gap-filling if it cannot be solved
+0267             failed=false;
+0268             try
+0269                 [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params);
+0270                 if exitFlag==-2
+0271                     EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n'];
+0272                     dispEM(EM,false);
+0273                 end
+0274             catch
+0275                 EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" could not be performed for any set of reactions\n'];
+0276                 dispEM(EM,false);
+0277                 failed=true;
+0278             end
+0279             if failed==false
+0280                 if ~isempty(newRxns)
+0281                     nAdded=nAdded+numel(newRxns);
+0282                     
+0283                     %Add the reactions to the base model. It is not correct
+0284                     %to use newModel directly, as it may contain
+0285                     %reactions/constraints that are specific to this task
+0286                     model=mergeModels({model,removeReactions(newModel,setdiff(newModel.rxns,newRxns),true,true)},'metNames',true);
+0287                     
+0288                     %Keep track of the added reactions
+0289                     addedRxns(ismember(refModel.rxns,newRxns),i)=true;
+0290                 end
+0291                 if printOutput==true
+0292                     fprintf(['[' taskStructure(i).id '] ' taskStructure(i).description ': Added ' num2str(numel(newRxns)) ' reaction(s), ' num2str(nAdded) ' reactions added in total\n']);
+0293                 end
+0294             end
+0295         else
+0296             if printOutput==true
+0297                 fprintf(['[' taskStructure(i).id '] ' taskStructure(i).description ': Added 0 reaction(s), ' num2str(nAdded) ' reactions added in total\n']);
+0298             end
+0299         end
+0300         supressWarnings=true;
+0301         
+0302         %Print the output if chosen
+0303         if taskStructure(i).printFluxes && printOutput
+0304             if ~isempty(sol.x)
+0305                 sol=solveLP(tModel,1);
+0306                 printFluxes(tModel,sol.x,false,10^-5,[],'%rxnID (%eqn):%flux\n');
+0307                 fprintf('\n');
+0308             else
+0309                 %If the problem wasn't solveable then the gap-filled model
+0310                 %should be used
+0311                 if failed==false
+0312                     sol=solveLP(newModel,1);
+0313                     printFluxes(newModel,sol.x,false,10^-5,[],'%rxnID (%eqn):%flux\n');
+0314                     fprintf('\n');
+0315                 end
+0316             end
+0317         end
+0318         
+0319         tModel=model;
+0320         %Since new mets are added by merging the new reactions and not only
+0321         %from the task sheet
+0322         modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
+0323     else
+0324         EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" is set as SHOULD FAIL. Such tasks cannot be modelled using this approach and the task is therefore ignored\n'];
+0325         dispEM(EM,false);
+0326     end
+0327 end
+0328 model.b(:,2) = [];  % resume field b
+0329 outModel=model;
+0330 end
    
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file
diff --git a/doc/core/gapReport.html b/doc/core/gapReport.html
index 84998030..3e8e78c5 100644
--- a/doc/core/gapReport.html
+++ b/doc/core/gapReport.html
@@ -71,7 +71,7 @@ 
    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
    
 This function calls:
 
+
  • canConsume	canConsume
  • canProduce	canProduce
  • checkProduction	checkProduction
  • fillGaps	fillGaps
  • getAllSubGraphs	getAllSubGraphs
  • haveFlux	haveFlux
  • removeReactions	removeReactions
  • simplifyModel	simplifyModel
    • 
 This function is called by:
 
@@ -180,7 +180,7 @@ 
    SOURCE CODE ^if ~isempty(templateModels)
 0100     fprintf('\n***Automated gap-filling\n');
-0101     [connectedFromTemplates, ~, addedFromTemplates]=fillGaps(model,templateModels);
+0101     [connectedFromTemplates, ~, addedFromTemplates]=fillGaps(model,templateModels);
 0102     t=templateModels{1}.id;
 0103     for i=2:numel(templateModels)
 0104         t=[t ', ' templateModels{i}.id];
diff --git a/doc/core/getExchangeRxns.html b/doc/core/getExchangeRxns.html
index 12648ae4..78d4e7a3 100644
--- a/doc/core/getExchangeRxns.html
+++ b/doc/core/getExchangeRxns.html
@@ -75,7 +75,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 </ul>
 This function is called by:
 <ul style=
-
  • consumeSomething	consumeSomething
  • fillGaps	fillGaps
  • makeSomething	makeSomething
  • predictLocalization	predictLocalization
  • removeBadRxns	removeBadRxns
  • runDynamicFBA	runDynamicFBA
  • setExchangeBounds	setExchangeBounds
    • 
+
  • consumeSomething	consumeSomething
  • fillGaps	fillGaps
  • makeSomething	makeSomething
  • predictLocalization	predictLocalization
  • removeBadRxns	removeBadRxns
  • runDynamicFBA	runDynamicFBA
  • setExchangeBounds	setExchangeBounds
    • 
 
 
 
diff --git a/doc/core/getMinNrFluxes.html b/doc/core/getMinNrFluxes.html
index 531c8df2..62a6847d 100644
--- a/doc/core/getMinNrFluxes.html
+++ b/doc/core/getMinNrFluxes.html
@@ -36,7 +36,7 @@ 
    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
 <li><a href=convertCharArray	convertCharArray
  • convertToIrrev	convertToIrrev
  • dispEM	dispEM
    • 
 This function is called by:
 
+
  • fillGaps	fillGaps
    • 
 
 
 
@@ -78,7 +78,7 @@ 
    SOURCE CODE ^%                 with the same number of elements as reactions in the model,
 0009 %                 of a vector of indexes for the reactions that should be
 0010 %                 minimized (optional, default model.rxns)
-0011 %   params        parameter structure as used by getMILPParams (optional)
+0011 %   params        *obsolete option*
 0012 %   scores        vector of weights for the reactions. Negative scores
 0013 %                 should not have flux. Positive scores are not possible in this
 0014 %                 implementation, and they are changed to max(scores(scores<0)).
diff --git a/doc/core/haveFlux.html b/doc/core/haveFlux.html
index 4fdbd8b6..4f10e9bc 100644
--- a/doc/core/haveFlux.html
+++ b/doc/core/haveFlux.html
@@ -56,7 +56,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 <li><a href=convertCharArray	convertCharArray
  • getIndexes	getIndexes
  • simplifyModel	simplifyModel
    • 
 This function is called by:
 
+
  • canConsume	canConsume
  • canProduce	canProduce
  • checkProduction	checkProduction
  • fillGaps	fillGaps
  • gapReport	gapReport
  • simplifyModel	simplifyModel
    • 
 
 
 
diff --git a/doc/core/makeSomething.html b/doc/core/makeSomething.html
index 16bf63ca..bbf9b56e 100644
--- a/doc/core/makeSomething.html
+++ b/doc/core/makeSomething.html
@@ -49,7 +49,7 @@ 
    DESCRIPTION ^SOURCE CODE ^%                   false)
 0022 %   allowExcretion  allow for excretion of all other metabolites (optional,
 0023 %                   default true)
-0024 %   params          parameter structure as used by getMILPParams (optional)
+0024 %   params          *obsolete option*
 0025 %   ignoreIntBounds    true if internal bounds (including reversibility)
 0026 %                   should be ignored. Exchange reactions are not affected.
 0027 %                   This can be used to find unbalanced solutions which are
diff --git a/doc/core/mergeModels.html b/doc/core/mergeModels.html
index 8264db25..b6d69f7c 100644
--- a/doc/core/mergeModels.html
+++ b/doc/core/mergeModels.html
@@ -53,7 +53,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 <li><a href=dispEM	dispEM
  • standardizeGrRules	standardizeGrRules
    • 
 This function is called by:
 
+
  • copyToComps	copyToComps
  • fillGaps	fillGaps
  • fitTasks	fitTasks
  • getModelFromHomology	getModelFromHomology
    • 
 
 
 
diff --git a/doc/core/parseTaskList.html b/doc/core/parseTaskList.html
index 71c1e0a9..091bf07f 100644
--- a/doc/core/parseTaskList.html
+++ b/doc/core/parseTaskList.html
@@ -144,7 +144,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 <li><a href=dispEM	dispEM
 This function is called by:
 
+
  • checkTasks	checkTasks
  • compareMultipleModels	compareMultipleModels
  • fitTasks	fitTasks
    • 
 
 
 
diff --git a/doc/core/printFluxes.html b/doc/core/printFluxes.html
index e9668d01..58765a14 100644
--- a/doc/core/printFluxes.html
+++ b/doc/core/printFluxes.html
@@ -70,7 +70,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 <li><a href=constructEquations	constructEquations
  • convertCharArray	convertCharArray
  • dispEM	dispEM
  • getElementalBalance	getElementalBalance
  • removeReactions	removeReactions
    • 
 This function is called by:
 
+
  • checkTasks	checkTasks
  • fitTasks	fitTasks
  • printModel	printModel
    • 
 
 
 
diff --git a/doc/core/randomSampling.html b/doc/core/randomSampling.html
index e089906f..17bf0e43 100644
--- a/doc/core/randomSampling.html
+++ b/doc/core/randomSampling.html
@@ -70,13 +70,18 @@ 
    DESCRIPTION ^SOURCE CODE ^% Output:
 0043 %   solutions               matrix with the solutions
 0044 %   goodRxns                double vector of indexes of those reactions
-0045 %                           that are not involved in loops and can be used
-0046 %                           as random objective functions
-0047 %
-0048 % The solutions are generated by maximizing (with random weights) for a
-0049 % random set of three reactions. For reversible reactions it randomly
-0050 % chooses between maximizing and minimizing.
-0051 %
-0052 % Usage: solutions = randomSampling(model, nSamples, replaceBoundsWithInf,...
-0053 %                       supressErrors, runParallel, goodRxns, minFlux)
-0054 
-0055 if nargin<2 | isempty(nSamples)
-0056     nSamples=1000;
-0057 end
-0058 if nargin<3 | isempty(replaceBoundsWithInf)
-0059     replaceBoundsWithInf=true;
-0060 end
-0061 if nargin<4 | isempty(supressErrors)
-0062     supressErrors=false;
-0063 end
-0064 if nargin<5 | isempty(runParallel)
-0065     runParallel=true;
-0066 end
-0067 if nargin<7 | isempty(minFlux)
-0068     minFlux=false;
-0069 end
-0070 
-0071 [ps, oldPoolAutoCreateSetting] = parallelPoolRAVEN(runParallel);
-0072 
-0073 nObjRxns=2; %Number of reactions in the objective function in each iteration
-0074 
-0075 %Simplify the model to speed stuff up a little. Keep original mapping
-0076 originalRxns=model.rxns;
-0077 model=simplifyModel(model,false,false,true,true);
-0078 
-0079 %Then change the bounds to +/- Inf. This is needed in order to not have
-0080 %loops in the solutions
-0081 if replaceBoundsWithInf==true
-0082     model.ub(model.ub==max(model.ub))=Inf;
-0083     if min(model.lb)<0 % Only negative lower bounds should be set to -Inf
-0084         model.lb(model.lb==min(model.lb))=-Inf;
-0085     end
-0086 end
-0087 
-0088 %Check that the model is feasible given the constraints
-0089 [sol,~]=solveLP(model);
-0090 if isempty(sol.x)
-0091     EM='The model has no feasible solution, likely due to incompatible constraints';
-0092     dispEM(EM);
-0093 elseif sol.f==0 && showProgress
-0094     warning('The model objective function cannot reach a non-zero value. This might be intended, so randomSampling will continue, but this could indicate problems with your model')
-0095 end
-0096 
-0097 [~,hsSol]=solveLP(model);
-0098 nRxns = numel(model.rxns);
-0099 %Reactions which can be involved in loops should not be optimized for.
-0100 %Check which reactions reach an arbitary high upper bound
-0101 if nargin<6 || isempty(goodRxns)
-0102     goodRxns = true(numel(model.rxns),numel(model.rxns));
-0103     goodRxns = num2cell(goodRxns,1);
-0104     PB = ProgressBar2(nRxns,'Prepare goodRxns not involved in loops','cli');
-0105     parfor (i=1:nRxns)
-0106         testModel=setParam(model,'eq',model.rxns(i),1000);
-0107         sol=solveLP(testModel,0,[],hsSol);
-0108         if ~isempty(sol.f)
-0109             notGood = abs(sol.x)>999;
-0110             goodRxns{i}(notGood)=false;
-0111         else
-0112             %If the reaction is reversible, also check in that direction
-0113             if model.rev(i)
-0114                 testModel=setParam(model,'eq',model.rxns(i),-1000);
-0115                 sol=solveLP(testModel,0,[],hsSol);
-0116                 if ~isempty(sol.f)
-0117                     notGood = abs(sol.x)>999;
-0118                     goodRxns{i}(notGood)=false;
-0119                 end
-0120             end
-0121         end
-0122         count(PB);
-0123     end
-0124     goodRxns = cell2mat(goodRxns);
-0125     goodRxns = find(prod(goodRxns,2));
-0126 end
-0127 
-0128 %Reserve space for a solution matrix
-0129 sols = zeros(numel(model.rxns),nSamples);
-0130 sols = num2cell(sols,1);
-0131 
-0132 %Main loop
-0133 if nSamples > 0
-0134 
-0135     PB = ProgressBar2(nSamples,'Performing random sampling','cli');
-0136     parfor i=1:nSamples
-0137         badSolutions = 1;
-0138         tmpModel = model;
-0139         while lt(0,badSolutions)
-0140             rxns = randsample(numel(goodRxns),nObjRxns);
-0141             tmpModel.c = zeros(numel(tmpModel.rxns),1);
-0142             multipliers = randsample([-1 1],nObjRxns,true);
-0143             multipliers(tmpModel.rev(goodRxns(rxns))==0) = 1;
-0144             tmpModel.c(goodRxns(rxns)) = rand(nObjRxns,1).*multipliers;
-0145             if true(minFlux)
-0146                 sol=solveLP(tmpModel,1,[],hsSol);
-0147             else
-0148                 sol=solveLP(tmpModel,0,[],hsSol);
-0149             end
-0150             if any(sol.x) && abs(sol.f)>10^-8
-0151                 sols{i} = sol.x;
-0152                 badSolutions = 0;
-0153             else
-0154                 badSolutions = badSolutions+1;
-0155                 %If it only finds bad solutions then throw an error.
-0156                 if badSolutions == 100 && supressErrors==false
-0157                     error('The program is having problems finding non-zero solutions (ignoring reactions that might be involved in loops). Review the constraints on your model. Set supressErrors to true to ignore this error');
-0158                 end
-0159             end
-0160         end
-0161         count(PB);
-0162     end
-0163 
-0164     %Map to original model
-0165     sols = cell2mat(sols);
-0166     [~, I]=ismember(model.rxns,originalRxns);
-0167     solutions=zeros(numel(originalRxns),nSamples);
-0168     solutions(I,:)=sols;
-0169 else
-0170     solutions=[];
-0171 end
-0172 % Reset original Parallel setting
-0173 ps.Pool.AutoCreate = oldPoolAutoCreateSetting;
-0174 end
+0045 %                           that are not involved in loops or always carry
+0046 %                           zero flux and can be used as random objective
+0047 %                           functions
+0048 %
+0049 % Note: The solutions are generated by maximizing (with random weights) for
+0050 % a random set of three reactions. For reversible reactions it randomly
+0051 % chooses between maximizing and minimizing.
+0052 %
+0053 % If the model is a GECKO v3+ ecModel, then usage_prot reactions are not
+0054 % selected for sampling, instead focusing on sampling from the metabolic
+0055 % aspects that form the solution space.
+0056 %
+0057 % Usage: solutions = randomSampling(model, nSamples, replaceBoundsWithInf,...
+0058 %                       supressErrors, runParallel, goodRxns, minFlux)
+0059 
+0060 if nargin<2 | isempty(nSamples)
+0061     nSamples=1000;
+0062 end
+0063 if nargin<3 | isempty(replaceBoundsWithInf)
+0064     replaceBoundsWithInf=true;
+0065 end
+0066 if nargin<4 | isempty(supressErrors)
+0067     supressErrors=false;
+0068 end
+0069 if nargin<5 | isempty(runParallel)
+0070     runParallel=true;
+0071 end
+0072 if nargin<7 | isempty(minFlux)
+0073     minFlux=false;
+0074 end
+0075 
+0076 [ps, oldPoolAutoCreateSetting] = parallelPoolRAVEN(runParallel);
+0077 
+0078 nObjRxns=2; %Number of reactions in the objective function in each iteration
+0079 
+0080 %Simplify the model to speed stuff up a little. Keep original mapping
+0081 originalRxns=model.rxns;
+0082 model=simplifyModel(model,false,false,true,true);
+0083 
+0084 %Then change the bounds to +/- Inf. This is needed in order to not have
+0085 %loops in the solutions
+0086 if replaceBoundsWithInf==true
+0087     model.ub(model.ub==max(model.ub))=Inf;
+0088     if min(model.lb)<0 % Only negative lower bounds should be set to -Inf
+0089         model.lb(model.lb==min(model.lb))=-Inf;
+0090     end
+0091 end
+0092 
+0093 %Check that the model is feasible given the constraints
+0094 [sol,~]=solveLP(model);
+0095 if isempty(sol.x)
+0096     EM='The model has no feasible solution, likely due to incompatible constraints';
+0097     dispEM(EM);
+0098 elseif sol.f==0 && showProgress
+0099     warning('The model objective function cannot reach a non-zero value. This might be intended, so randomSampling will continue, but this could indicate problems with your model')
+0100 end
+0101 
+0102 nRxns = numel(model.rxns);
+0103 %Reactions which can be involved in loops should not be optimized for.
+0104 %Check which reactions reach an arbitary high upper bound
+0105 if nargin<6 || isempty(goodRxns)
+0106     revRxns = transpose(find(model.lb < 0));
+0107     fwdRxns = transpose(find(model.ub > 0));
+0108     rxnList = [fwdRxns,revRxns];
+0109     objList = [ones(numel(fwdRxns),1);-ones(numel(revRxns),1)];
+0110     testSol = zeros(numel(objList),1);
+0111     PB = ProgressBar2(nRxns,'Prepare goodRxns not involved in loops','cli');
+0112 
+0113     parfor i = 1:numel(objList)
+0114         testModel=setParam(model,'obj',rxnList(i),objList(i));
+0115         sol=solveLP(testModel,0);
+0116         if ~isempty(sol.f)
+0117             testSol(i) = sol.x(rxnList(i));
+0118         end
+0119         count(PB);
+0120     end
+0121     testSol  = testSol.*objList;
+0122     goodRxns = true(nRxns,1);
+0123     % Filter out reactions that can reach (-)1000 (= involved in loop)
+0124     goodRxns(fwdRxns(testSol(1:numel(fwdRxns)) > 999)) = false;
+0125     goodRxns(revRxns(testSol(numel(fwdRxns)+1:end) > 999)) = false;
+0126     testSol(revRxns) = testSol(revRxns) + testSol(numel(fwdRxns)+1:end);
+0127     % Filter out reactions that cannot carry flux
+0128     goodRxns(testSol(1:nRxns) == 0) = false;
+0129     % In ecModels do not sample from usage_prot reactions
+0130     if isfield(model,'ec')
+0131         usageRxns = startsWith(model.rxns,'usage_prot_');
+0132         goodRxns(usageRxns) = false;
+0133     end
+0134     goodRxns = find(goodRxns);
+0135 end
+0136 
+0137 %Reserve space for a solution matrix
+0138 sols = zeros(numel(model.rxns),nSamples);
+0139 sols = num2cell(sols,1);
+0140 
+0141 %Main loop
+0142 if nSamples > 0
+0143 
+0144     PB = ProgressBar2(nSamples,'Performing random sampling','cli');
+0145     parfor i=1:nSamples
+0146         badSolutions = 1;
+0147         tmpModel = model;
+0148         while lt(0,badSolutions)
+0149             rxns = randsample(numel(goodRxns),nObjRxns);
+0150             tmpModel.c = zeros(numel(tmpModel.rxns),1);
+0151             
+0152             multipliers = randsample([-1 1],nObjRxns,true);
+0153             multipliers(tmpModel.ub(goodRxns(rxns))<0) = 1;
+0154             
+0155             tmpModel.c(goodRxns(rxns)) = rand(nObjRxns,1).*multipliers;
+0156 
+0157             if true(minFlux)
+0158                 sol=solveLP(tmpModel,1);
+0159             else
+0160                 sol=solveLP(tmpModel,0);
+0161             end
+0162             if any(sol.x) && abs(sol.f)>10^-8
+0163                 sols{i} = sol.x;
+0164                 badSolutions = 0;
+0165             else
+0166                 badSolutions = badSolutions+1;
+0167                 %If it only finds bad solutions then throw an error.
+0168                 if badSolutions == 100 && supressErrors==false
+0169                     error('The program is having problems finding non-zero solutions (ignoring reactions that might be involved in loops). Review the constraints on your model. Set supressErrors to true to ignore this error');
+0170                 end
+0171             end
+0172         end
+0173         count(PB);
+0174     end
 0175 
-0176 %To use instead of the normal Matlab randsample function. This is in order
-0177 %to not depend on the Matlab statistical toolbox.
-0178 function I=randsample(n,k,replacement)
-0179 if nargin<3
-0180     replacement=false;
-0181 end
-0182 %n can be a integer, which leads to I being sampled from 1:n, or it can be
-0183 %a population to sample from.
-0184 if numel(n)==1 && isnumeric(n)
-0185     n=1:n;
+0176     %Map to original model
+0177     sols = cell2mat(sols);
+0178     [~, I]=ismember(model.rxns,originalRxns);
+0179     solutions=zeros(numel(originalRxns),nSamples);
+0180     solutions(I,:)=sols;
+0181 else
+0182     solutions=[];
+0183 end
+0184 % Reset original Parallel setting
+0185 ps.Pool.AutoCreate = oldPoolAutoCreateSetting;
 0186 end
-0187 %Loop and get random numbers until the list is unique. This is only a good
-0188 %option is the number of samples is small compared to the population. There
-0189 %are several checks that should be made here, for example regarding size
-0190 %and that the number of samples is <=population size if replacement==false.
-0191 %This is not the case in randomSampling, so such checks are ignored
-0192 while true
-0193     J=randi(numel(n),[k,1]);
-0194     if replacement==true || numel(J)==numel(unique(J))
-0195         I=n(J);
-0196         break;
-0197     end
+0187 
+0188 %To use instead of the normal Matlab randsample function. This is in order
+0189 %to not depend on the Matlab statistical toolbox.
+0190 function I=randsample(n,k,replacement)
+0191 if nargin<3
+0192     replacement=false;
+0193 end
+0194 %n can be a integer, which leads to I being sampled from 1:n, or it can be
+0195 %a population to sample from.
+0196 if numel(n)==1 && isnumeric(n)
+0197     n=1:n;
 0198 end
-0199 I=I(:);
-0200 end
    +0199 %Loop and get random numbers until the list is unique. This is only a good +0200 %option is the number of samples is small compared to the population. There +0201 %are several checks that should be made here, for example regarding size +0202 %and that the number of samples is <=population size if replacement==false. +0203 %This is not the case in randomSampling, so such checks are ignored +0204 while true +0205 J=randi(numel(n),[k,1]); +0206 if replacement==true || numel(J)==numel(unique(J)) +0207 I=n(J); +0208 break; +0209 end +0210 end +0211 I=I(:); +0212 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/core/removeReactions.html b/doc/core/removeReactions.html index efc52428..0f658fe0 100644 --- a/doc/core/removeReactions.html +++ b/doc/core/removeReactions.html @@ -55,7 +55,7 @@

    CROSS-REFERENCE INFORMATION ^ <li><a href=convertCharArray convertCharArray
  • getIndexes getIndexes
  • removeMets removeMets
    • This function is called by: +
  • changeRxns changeRxns
  • contractModel contractModel
  • copyToComps copyToComps
  • fillGaps fillGaps
  • findGeneDeletions findGeneDeletions
  • fitTasks fitTasks
  • gapReport gapReport
  • getElementalBalance getElementalBalance
  • getModelFromHomology getModelFromHomology
  • mergeCompartments mergeCompartments
  • predictLocalization predictLocalization
  • printFluxes printFluxes
  • removeBadRxns removeBadRxns
  • removeGenes removeGenes
  • removeMets removeMets
  • simplifyModel simplifyModel
    • diff --git a/doc/core/setParam.html b/doc/core/setParam.html index 2aee63e2..73b63b7e 100644 --- a/doc/core/setParam.html +++ b/doc/core/setParam.html @@ -64,7 +64,7 @@

    CROSS-REFERENCE INFORMATION ^ <li><a href=convertCharArray convertCharArray
  • dispEM dispEM
    • This function is called by: +
  • FSEOF FSEOF
  • checkRxn checkRxn
  • checkTasks checkTasks
  • consumeSomething consumeSomething
  • fitTasks fitTasks
  • getEssentialRxns getEssentialRxns
  • makeSomething makeSomething
  • randomSampling randomSampling
  • runPhenotypePhasePlane runPhenotypePhasePlane
  • runProductionEnvelope runProductionEnvelope
  • runRobustnessAnalysis runRobustnessAnalysis
  • runSimpleOptKnock runSimpleOptKnock
  • setExchangeBounds setExchangeBounds
    • @@ -195,7 +195,12 @@

    SOURCE CODE ^end 0124 end 0125 end -0126 end +0126 if any(ismember(paramType,{'lb','ub','unc'})) +0127 invalidBound = model.lb(indexes) > model.ub(indexes); +0128 if any(invalidBound) +0129 error(['Invalid set of bounds for reaction(s): ', strjoin(model.rxns(indexes(invalidBound)),', '), '.']) +0130 end +0131 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/core/simplifyModel.html b/doc/core/simplifyModel.html index f6f4baf6..b47bde2d 100644 --- a/doc/core/simplifyModel.html +++ b/doc/core/simplifyModel.html @@ -70,7 +70,7 @@

    CROSS-REFERENCE INFORMATION ^ <li><a href=contractModel contractModel
  • convertCharArray convertCharArray
  • convertToIrrev convertToIrrev
  • dispEM dispEM
  • getAllowedBounds getAllowedBounds
  • getIndexes getIndexes
  • haveFlux haveFlux
  • removeMets removeMets
  • removeReactions removeReactions
    • This function is called by: +
  • fillGaps fillGaps
  • findGeneDeletions findGeneDeletions
  • gapReport gapReport
  • haveFlux haveFlux
  • printModelStats printModelStats
  • randomSampling randomSampling
    • diff --git a/doc/io/importModel.html b/doc/io/importModel.html index 05f75018..38d2492f 100644 --- a/doc/io/importModel.html +++ b/doc/io/importModel.html @@ -265,7 +265,7 @@

    SOURCE CODE ^if isfield(modelSBML.compartment,'sboTerm') && numel(unique([modelSBML.compartment.sboTerm])) == 1 +0158 if isfield(modelSBML.compartment,'sboTerm') && isscalar(unique([modelSBML.compartment.sboTerm])) 0159 %If all the SBO terms are identical, don't add them to compMiriams 0160 modelSBML.compartment = rmfield(modelSBML.compartment,'sboTerm'); 0161 end @@ -401,7 +401,7 @@

    SOURCE CODE ^... 0292 formula(compositionIndexes(1)+1:compositionIndexes(2)-1); 0293 else -0294 if numel(compositionIndexes)==1 +0294 if isscalar(compositionIndexes) 0295 %Probably a simple molecule which can have only 0296 %one conformation 0297 metaboliteFormula{numel(metaboliteFormula)+1,1}=... @@ -532,7 +532,7 @@

    SOURCE CODE ^end 0424 -0425 if isfield(modelSBML.reaction,'sboTerm') && numel(unique([modelSBML.reaction.sboTerm])) == 1 +0425 if isfield(modelSBML.reaction,'sboTerm') && isscalar(unique([modelSBML.reaction.sboTerm])) 0426 %If all the SBO terms are identical, don't add them to rxnMiriams 0427 modelSBML.reaction = rmfield(modelSBML.reaction,'sboTerm'); 0428 end @@ -543,7 +543,7 @@

    SOURCE CODE ^%so, then jump to the next reaction. This is because I get the genes 0434 %for complexes from the names and not from the reactions that create 0435 %them. This only applies to the non-COBRA format -0436 if numel(modelSBML.reaction(i).product)==1 +0436 if isscalar(modelSBML.reaction(i).product) 0437 if length(modelSBML.reaction(i).product(1).species)>=3 0438 if strcmp(modelSBML.reaction(i).product(1).species(1:3),'Cx_')==true 0439 continue; @@ -565,701 +565,695 @@

    SOURCE CODE ^if isfield(modelSBML.reaction(i),'fbc_lowerFluxBound') 0456 lb=modelSBML.reaction(i).fbc_lowerFluxBound; 0457 ub=modelSBML.reaction(i).fbc_upperFluxBound; -0458 for n=1:numel(parameter.value) -0459 lb=regexprep(lb,parameter.name(n),num2str(parameter.value{n})); -0460 ub=regexprep(ub,parameter.name(n),num2str(parameter.value{n})); -0461 end -0462 if isempty(lb) -0463 lb='-Inf'; -0464 end -0465 if isempty(ub) -0466 ub='Inf'; -0467 end -0468 reactionLB(counter)=str2num(lb); -0469 reactionUB(counter)=str2num(ub); -0470 %The order of these parameters should not be hard coded -0471 elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter') -0472 reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value; -0473 reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value; -0474 reactionObjective(counter)=modelSBML.reaction(i).kineticLaw.parameter(3).value; -0475 else -0476 if reactionReversibility(counter)==true -0477 reactionLB(counter)=-inf; -0478 else -0479 reactionLB(counter)=0; +0458 [~,fluxBoundIdx] = ismember({lb,ub},parameter.name); +0459 reactionLB(counter) = parameter.value{fluxBoundIdx(1)}; +0460 reactionUB(counter) = parameter.value{fluxBoundIdx(2)}; +0461 elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter') +0462 reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value; +0463 reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value; +0464 reactionObjective(counter)=modelSBML.reaction(i).kineticLaw.parameter(3).value; +0465 else +0466 if reactionReversibility(counter)==true +0467 reactionLB(counter)=-inf; +0468 else +0469 reactionLB(counter)=0; +0470 end +0471 reactionUB(counter)=inf; +0472 reactionObjective(counter)=0; +0473 end +0474 +0475 %Find the associated gene if available +0476 %If FBC, get gene association data from corresponding fields +0477 if isfield(modelSBML.reaction(i),'fbc_geneProductAssociation') +0478 if ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation) && ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association) +0479 grRules{counter}=modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association; 0480 end -0481 reactionUB(counter)=inf; -0482 reactionObjective(counter)=0; -0483 end -0484 -0485 %Find the associated gene if available -0486 %If FBC, get gene association data from corresponding fields -0487 if isfield(modelSBML.reaction(i),'fbc_geneProductAssociation') -0488 if ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation) && ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association) -0489 grRules{counter}=modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association; -0490 end -0491 elseif isfield(modelSBML.reaction(i),'notes') -0492 %This section was previously executed only if isSBML2COBRA is true. Now -0493 %it will be executed, if 'GENE_ASSOCIATION' is found in -0494 %modelSBML.reaction(i).notes -0495 if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION') -0496 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION'); -0497 elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION') -0498 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION'); -0499 else -0500 geneAssociation=''; -0501 end -0502 if ~isempty(geneAssociation) -0503 %This adds the grRules. The gene list and rxnGeneMat are created -0504 %later -0505 grRules{counter}=geneAssociation; -0506 end -0507 end -0508 if isempty(grRules{counter}) && ~isempty(modelSBML.reaction(i).modifier) -0509 rules=''; -0510 for j=1:numel(modelSBML.reaction(i).modifier) -0511 modifier=modelSBML.reaction(i).modifier(j).species; -0512 if ~isempty(modifier) -0513 if strcmpi(modifier(1:2),'E_') -0514 index=find(strcmp(modifier,geneIDs)); -0515 %This should be unique and in the geneIDs list, -0516 %otherwise something is wrong -0517 if numel(index)~=1 -0518 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; -0519 dispEM(EM); -0520 end -0521 if ~isempty(rules) -0522 rules=[rules ' or (' geneNames{index} ')']; -0523 else -0524 rules=['(' geneNames{index} ')']; -0525 end -0526 elseif strcmp(modifier(1:2),'s_') -0527 index=find(strcmp(modifier,metaboliteIDs)); -0528 %This should be unique and in the geneIDs list, -0529 %otherwise something is wrong -0530 if numel(index)~=1 -0531 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; -0532 dispEM(EM); -0533 end -0534 if ~isempty(rules) -0535 rules=[rules ' or (' metaboliteIDs{index} ')']; -0536 else -0537 rules=['(' metaboliteIDs{index} ')']; -0538 end -0539 else -0540 %It seems to be a complex. Add the corresponding -0541 %genes from the name of the complex (not the -0542 %reaction that creates it) -0543 index=find(strcmp(modifier,complexIDs)); -0544 if numel(index)==1 -0545 if ~isempty(rules) -0546 rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')']; -0547 else -0548 rules=['(' strrep(complexNames{index},':',' and ') ')']; -0549 end -0550 else -0551 %Could not find a complex -0552 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; -0553 dispEM(EM); -0554 end -0555 end -0556 end -0557 end -0558 grRules{counter}=rules; -0559 grRulesFromModifier{counter}=rules;%Backup copy for grRules, useful to parse Yeast 7.6 -0560 end -0561 -0562 %Add reaction compartment -0563 if isfield(modelSBML.reaction(i),'compartment') -0564 if ~isempty(modelSBML.reaction(i).compartment) -0565 rxnComp=modelSBML.reaction(i).compartment; -0566 else -0567 rxnComp=''; -0568 end -0569 elseif isfield(modelSBML.reaction(i),'notes') -0570 rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT'); -0571 end -0572 if ~isempty(rxnComp) -0573 %Find it in the compartment list -0574 [~, J]=ismember(rxnComp,compartmentIDs); -0575 rxnComps(counter)=J; -0576 end -0577 -0578 -0579 miriamStruct=parseMiriam(modelSBML.reaction(i).annotation); -0580 rxnMiriams{counter}=miriamStruct; -0581 if isfield(modelSBML.reaction(i),'notes') -0582 subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM')); -0583 subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = []; -0584 if strfind(modelSBML.reaction(i).notes,'Confidence Level') -0585 confScore = parseNote(modelSBML.reaction(i).notes,'Confidence Level'); -0586 if isempty(confScore) -0587 confScore = 0; -0588 end -0589 rxnconfidencescores(counter)=str2double(confScore); -0590 end -0591 rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS'); -0592 rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES'); -0593 end -0594 -0595 %Get SBO terms -0596 if isfield(modelSBML.reaction(i),'sboTerm') && ~(modelSBML.reaction(i).sboTerm==-1) -0597 rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm); -0598 end -0599 -0600 %Get ec-codes -0601 eccode=''; -0602 if ~isempty(modelSBML.reaction(i).annotation) -0603 if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code') -0604 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code'); -0605 elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code') -0606 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code'); -0607 elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code') -0608 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code'); -0609 end -0610 elseif isfield(modelSBML.reaction(i),'notes') -0611 if strfind(modelSBML.reaction(i).notes,'EC Number') -0612 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')]; -0613 elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS') -0614 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')]; -0615 end -0616 end -0617 eccodes{counter}=eccode; -0618 -0619 %Add all reactants -0620 for j=1:numel(modelSBML.reaction(i).reactant) -0621 %Get the index of the metabolite in metaboliteIDs. External -0622 %metabolites will be removed at a later stage -0623 metIndex=find(strcmp(modelSBML.reaction(i).reactant(j).species,metaboliteIDs),1); -0624 if isempty(metIndex) -0625 EM=['Could not find metabolite ' modelSBML.reaction(i).reactant(j).species ' in reaction ' reactionIDs{counter}]; -0626 dispEM(EM); -0627 end -0628 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).reactant(j).stoichiometry*-1; -0629 end -0630 -0631 %Add all products -0632 for j=1:numel(modelSBML.reaction(i).product) -0633 %Get the index of the metabolite in metaboliteIDs. -0634 metIndex=find(strcmp(modelSBML.reaction(i).product(j).species,metaboliteIDs),1); -0635 if isempty(metIndex) -0636 EM=['Could not find metabolite ' modelSBML.reaction(i).product(j).species ' in reaction ' reactionIDs{counter}]; -0637 dispEM(EM); -0638 end -0639 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).product(j).stoichiometry; -0640 end -0641 end -0642 -0643 %if FBC, objective function is separately defined. Multiple objective -0644 %functions can be defined, one is set as active -0645 if isfield(modelSBML, 'fbc_activeObjective') -0646 obj=modelSBML.fbc_activeObjective; -0647 for i=1:numel(modelSBML.fbc_objective) -0648 if strcmp(obj,modelSBML.fbc_objective(i).fbc_id) -0649 if ~isempty(modelSBML.fbc_objective(i).fbc_fluxObjective) -0650 rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction; -0651 idx=find(ismember(reactionIDs,rxn)); -0652 reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient; -0653 end -0654 end -0655 end -0656 end -0657 -0658 %subSystems can be stored as groups instead of in annotations -0659 if isfield(modelSBML,'groups_group') -0660 for i=1:numel(modelSBML.groups_group) -0661 groupreactions={modelSBML.groups_group(i).groups_member(:).groups_idRef}; -0662 [~, idx] = ismember(groupreactions, reactionIDs); -0663 if any(idx) -0664 for j=1:numel(idx) -0665 if isempty(subsystems{idx(j)}) % First subsystem -0666 subsystems{idx(j)} = {modelSBML.groups_group(i).groups_name}; -0667 else % Consecutive subsystems: concatenate -0668 subsystems{idx(j)} = horzcat(subsystems{idx(j)}, modelSBML.groups_group(i).groups_name); -0669 end -0670 end -0671 end -0672 end -0673 end -0674 -0675 %Shrink the structures if complex-forming reactions had to be skipped -0676 reactionNames=reactionNames(1:counter); -0677 reactionIDs=reactionIDs(1:counter); -0678 subsystems=subsystems(1:counter); -0679 eccodes=eccodes(1:counter); -0680 rxnconfidencescores=rxnconfidencescores(1:counter); -0681 rxnreferences=rxnreferences(1:counter); -0682 rxnnotes=rxnnotes(1:counter); -0683 grRules=grRules(1:counter); -0684 rxnMiriams=rxnMiriams(1:counter); -0685 reactionReversibility=reactionReversibility(1:counter); -0686 reactionUB=reactionUB(1:counter); -0687 reactionLB=reactionLB(1:counter); -0688 reactionObjective=reactionObjective(1:counter); -0689 S=S(:,1:counter); -0690 -0691 model.name=modelSBML.name; -0692 model.id=modelSBML.id; -0693 model.rxns=reactionIDs; -0694 model.mets=metaboliteIDs; -0695 model.S=sparse(S); -0696 model.lb=reactionLB; -0697 model.ub=reactionUB; -0698 model.rev=reactionReversibility; -0699 model.c=reactionObjective; -0700 model.b=zeros(numel(metaboliteIDs),1); -0701 model.comps=compartmentIDs; -0702 model.compNames=compartmentNames; -0703 model.rxnConfidenceScores=rxnconfidencescores; -0704 model.rxnReferences=rxnreferences; -0705 model.rxnNotes=rxnnotes; -0706 -0707 %Load annotation if available. If there are several authors, only the first -0708 %author credentials are imported -0709 if isfield(modelSBML,'annotation') -0710 endString='</'; -0711 I=strfind(modelSBML.annotation,endString); -0712 J=strfind(modelSBML.annotation,'<vCard:Family>'); -0713 if any(J) -0714 model.annotation.familyName=modelSBML.annotation(J(1)+14:I(find(I>J(1),1))-1); -0715 end -0716 J=strfind(modelSBML.annotation,'<vCard:Given>'); -0717 if any(J) -0718 model.annotation.givenName=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); -0719 end -0720 J=strfind(modelSBML.annotation,'<vCard:EMAIL>'); -0721 if any(J) -0722 model.annotation.email=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); -0723 end -0724 J=strfind(modelSBML.annotation,'<vCard:Orgname>'); -0725 if any(J) -0726 model.annotation.organization=modelSBML.annotation(J(1)+15:I(find(I>J(1),1))-1); -0727 end -0728 endString='"/>'; -0729 I=strfind(modelSBML.annotation,endString); -0730 if strfind(modelSBML.annotation,'"urn:miriam:') -0731 J=strfind(modelSBML.annotation,'"urn:miriam:'); -0732 if any(J) -0733 model.annotation.taxonomy=modelSBML.annotation(J+12:I(find(I>J,1))-1); -0734 end -0735 else -0736 J=strfind(modelSBML.annotation,'"http://identifiers.org/'); -0737 if any(J) -0738 model.annotation.taxonomy=modelSBML.annotation(J+24:I(find(I>J,1))-1); -0739 else -0740 J=strfind(modelSBML.annotation,'"https://identifiers.org/'); -0741 if any(J) -0742 model.annotation.taxonomy=modelSBML.annotation(J+25:I(find(I>J,1))-1); -0743 end -0744 end -0745 end -0746 end -0747 if isfield(modelSBML,'notes') -0748 startString=strfind(modelSBML.notes,'xhtml">'); -0749 endString=strfind(modelSBML.notes,'</body>'); -0750 if any(startString) && any(endString) -0751 model.annotation.note=modelSBML.notes(startString+7:endString-1); -0752 model.annotation.note=regexprep(model.annotation.note,'<p>|</p>',''); -0753 model.annotation.note=strtrim(model.annotation.note); -0754 if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') -0755 model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel -0756 end -0757 end -0758 end -0759 -0760 if any(~cellfun(@isempty,compartmentOutside)) -0761 model.compOutside=compartmentOutside; -0762 end -0763 -0764 model.rxnNames=reactionNames; -0765 model.metNames=metaboliteNames; -0766 -0767 %Match the compartments for metabolites -0768 [~, J]=ismember(metaboliteCompartments,model.comps); -0769 model.metComps=J; -0770 -0771 %If any genes have been loaded (only for the new format) -0772 if ~isempty(geneNames) -0773 %In some rare cases geneNames may not necessarily be used in grRules. -0774 %That is true for Yeast 7.6. It's therefore important to change gene -0775 %systematic names to geneIDs in sophisticated way. Gene systematic -0776 %names are not unique, since exactly the same name may be in different -0777 %compartments -0778 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) -0779 geneShortNames=geneNames; -0780 %geneShortNames contain compartments as well, so these are removed -0781 geneShortNames=regexprep(geneShortNames,' \[.+$',''); -0782 %grRules obtained from modifier fields contain geneNames. These are -0783 %changed into geneIDs. grRulesFromModifier is a good way to have -0784 %geneIDs and rxns association when it's important to resolve -0785 %systematic name ambiguities -0786 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0787 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0788 -0789 %Yeast 7.6 contains several metabolites, which were used in gene -0790 %associations. For that reason, the list of species ID is created -0791 %and we then check whether any of them have kegg.genes annotation -0792 %thereby obtaining systematic gene names -0793 geneShortNames=vertcat(geneShortNames,metaboliteNames); -0794 geneIDs=vertcat(geneIDs,metaboliteIDs); -0795 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); -0796 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); -0797 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); -0798 -0799 %Now we retain information for only these entries, which have -0800 %kegg.genes annotation -0801 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); -0802 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); -0803 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); -0804 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); -0805 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); -0806 %Now we reorder geneIDs and geneSystNames by geneSystNames string -0807 %length -0808 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later -0809 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); -0810 geneIDs = geneIDs(Indx); -0811 geneSystNames = geneSystNames(Indx); -0812 for i=1:numel(geneSystNames) -0813 for j=1:numel(grRules) -0814 if strfind(grRules{j},geneSystNames{i}) -0815 if ~isempty(grRules{j}) -0816 if sum(ismember(geneSystNames,geneSystNames{i}))==1 -0817 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); -0818 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 -0819 counter=0; -0820 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); -0821 for k=1:numel(ovrlpIDs) -0822 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) -0823 counter=counter+1; -0824 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); -0825 end -0826 if counter>1 -0827 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; -0828 dispEM(EM); -0829 end -0830 end -0831 end -0832 end -0833 end -0834 end -0835 end -0836 end -0837 model.genes=geneNames; -0838 model.grRules=grRules; -0839 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0840 model.grRules = grRules; -0841 model.rxnGeneMat = rxnGeneMat; -0842 -0843 %Match the compartments for genes -0844 [~, J]=ismember(geneCompartments,model.comps); -0845 model.geneComps=J; -0846 else -0847 if ~all(cellfun(@isempty,grRules)) -0848 %If fbc_geneProduct exists, follow the specified gene order, such -0849 %that matching geneShortNames in function below will work -0850 if isfield(modelSBML,'fbc_geneProduct') -0851 genes={modelSBML.fbc_geneProduct.fbc_id}; -0852 -0853 %Get gene Miriams if they were not retrieved above (this occurs -0854 %when genes are stored as fbc_geneProduct instead of species) -0855 if isempty(geneMiriams) -0856 geneMiriams = cell(numel(genes),1); -0857 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && numel(unique([modelSBML.fbc_geneProduct.sboTerm])) == 1 -0858 %If all the SBO terms are identical, don't add them to geneMiriams -0859 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); +0481 elseif isfield(modelSBML.reaction(i),'notes') +0482 %This section was previously executed only if isSBML2COBRA is true. Now +0483 %it will be executed, if 'GENE_ASSOCIATION' is found in +0484 %modelSBML.reaction(i).notes +0485 if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION') +0486 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION'); +0487 elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION') +0488 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION'); +0489 else +0490 geneAssociation=''; +0491 end +0492 if ~isempty(geneAssociation) +0493 %This adds the grRules. The gene list and rxnGeneMat are created +0494 %later +0495 grRules{counter}=geneAssociation; +0496 end +0497 end +0498 if isempty(grRules{counter}) && ~isempty(modelSBML.reaction(i).modifier) +0499 rules=''; +0500 for j=1:numel(modelSBML.reaction(i).modifier) +0501 modifier=modelSBML.reaction(i).modifier(j).species; +0502 if ~isempty(modifier) +0503 if strcmpi(modifier(1:2),'E_') +0504 index=find(strcmp(modifier,geneIDs)); +0505 %This should be unique and in the geneIDs list, +0506 %otherwise something is wrong +0507 if numel(index)~=1 +0508 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; +0509 dispEM(EM); +0510 end +0511 if ~isempty(rules) +0512 rules=[rules ' or (' geneNames{index} ')']; +0513 else +0514 rules=['(' geneNames{index} ')']; +0515 end +0516 elseif strcmp(modifier(1:2),'s_') +0517 index=find(strcmp(modifier,metaboliteIDs)); +0518 %This should be unique and in the geneIDs list, +0519 %otherwise something is wrong +0520 if numel(index)~=1 +0521 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; +0522 dispEM(EM); +0523 end +0524 if ~isempty(rules) +0525 rules=[rules ' or (' metaboliteIDs{index} ')']; +0526 else +0527 rules=['(' metaboliteIDs{index} ')']; +0528 end +0529 else +0530 %It seems to be a complex. Add the corresponding +0531 %genes from the name of the complex (not the +0532 %reaction that creates it) +0533 index=find(strcmp(modifier,complexIDs)); +0534 if isscalar(index) +0535 if ~isempty(rules) +0536 rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')']; +0537 else +0538 rules=['(' strrep(complexNames{index},':',' and ') ')']; +0539 end +0540 else +0541 %Could not find a complex +0542 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; +0543 dispEM(EM); +0544 end +0545 end +0546 end +0547 end +0548 grRules{counter}=rules; +0549 grRulesFromModifier{counter}=rules;%Backup copy for grRules, useful to parse Yeast 7.6 +0550 end +0551 +0552 %Add reaction compartment +0553 if isfield(modelSBML.reaction(i),'compartment') +0554 if ~isempty(modelSBML.reaction(i).compartment) +0555 rxnComp=modelSBML.reaction(i).compartment; +0556 else +0557 rxnComp=''; +0558 end +0559 elseif isfield(modelSBML.reaction(i),'notes') +0560 rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT'); +0561 end +0562 if ~isempty(rxnComp) +0563 %Find it in the compartment list +0564 [~, J]=ismember(rxnComp,compartmentIDs); +0565 rxnComps(counter)=J; +0566 end +0567 +0568 +0569 miriamStruct=parseMiriam(modelSBML.reaction(i).annotation); +0570 rxnMiriams{counter}=miriamStruct; +0571 if isfield(modelSBML.reaction(i),'notes') +0572 subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM')); +0573 subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = []; +0574 if strfind(modelSBML.reaction(i).notes,'Confidence Level') +0575 confScore = parseNote(modelSBML.reaction(i).notes,'Confidence Level'); +0576 if isempty(confScore) +0577 confScore = 0; +0578 end +0579 rxnconfidencescores(counter)=str2double(confScore); +0580 end +0581 rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS'); +0582 rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES'); +0583 end +0584 +0585 %Get SBO terms +0586 if isfield(modelSBML.reaction(i),'sboTerm') && ~(modelSBML.reaction(i).sboTerm==-1) +0587 rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm); +0588 end +0589 +0590 %Get ec-codes +0591 eccode=''; +0592 if ~isempty(modelSBML.reaction(i).annotation) +0593 if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code') +0594 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code'); +0595 elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code') +0596 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code'); +0597 elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code') +0598 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code'); +0599 end +0600 elseif isfield(modelSBML.reaction(i),'notes') +0601 if strfind(modelSBML.reaction(i).notes,'EC Number') +0602 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')]; +0603 elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS') +0604 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')]; +0605 end +0606 end +0607 eccodes{counter}=eccode; +0608 +0609 %Add all reactants +0610 for j=1:numel(modelSBML.reaction(i).reactant) +0611 %Get the index of the metabolite in metaboliteIDs. External +0612 %metabolites will be removed at a later stage +0613 metIndex=find(strcmp(modelSBML.reaction(i).reactant(j).species,metaboliteIDs),1); +0614 if isempty(metIndex) +0615 EM=['Could not find metabolite ' modelSBML.reaction(i).reactant(j).species ' in reaction ' reactionIDs{counter}]; +0616 dispEM(EM); +0617 end +0618 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).reactant(j).stoichiometry*-1; +0619 end +0620 +0621 %Add all products +0622 for j=1:numel(modelSBML.reaction(i).product) +0623 %Get the index of the metabolite in metaboliteIDs. +0624 metIndex=find(strcmp(modelSBML.reaction(i).product(j).species,metaboliteIDs),1); +0625 if isempty(metIndex) +0626 EM=['Could not find metabolite ' modelSBML.reaction(i).product(j).species ' in reaction ' reactionIDs{counter}]; +0627 dispEM(EM); +0628 end +0629 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).product(j).stoichiometry; +0630 end +0631 end +0632 +0633 %if FBC, objective function is separately defined. Multiple objective +0634 %functions can be defined, one is set as active +0635 if isfield(modelSBML, 'fbc_activeObjective') +0636 obj=modelSBML.fbc_activeObjective; +0637 for i=1:numel(modelSBML.fbc_objective) +0638 if strcmp(obj,modelSBML.fbc_objective(i).fbc_id) +0639 if ~isempty(modelSBML.fbc_objective(i).fbc_fluxObjective) +0640 rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction; +0641 idx=ismember(reactionIDs,rxn); +0642 reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient; +0643 end +0644 end +0645 end +0646 end +0647 +0648 %subSystems can be stored as groups instead of in annotations +0649 if isfield(modelSBML,'groups_group') +0650 for i=1:numel(modelSBML.groups_group) +0651 groupreactions={modelSBML.groups_group(i).groups_member(:).groups_idRef}; +0652 [~, idx] = ismember(groupreactions, reactionIDs); +0653 if any(idx) +0654 for j=1:numel(idx) +0655 if isempty(subsystems{idx(j)}) % First subsystem +0656 subsystems{idx(j)} = {modelSBML.groups_group(i).groups_name}; +0657 else % Consecutive subsystems: concatenate +0658 subsystems{idx(j)} = horzcat(subsystems{idx(j)}, modelSBML.groups_group(i).groups_name); +0659 end +0660 end +0661 end +0662 end +0663 end +0664 +0665 %Shrink the structures if complex-forming reactions had to be skipped +0666 reactionNames=reactionNames(1:counter); +0667 reactionIDs=reactionIDs(1:counter); +0668 subsystems=subsystems(1:counter); +0669 eccodes=eccodes(1:counter); +0670 rxnconfidencescores=rxnconfidencescores(1:counter); +0671 rxnreferences=rxnreferences(1:counter); +0672 rxnnotes=rxnnotes(1:counter); +0673 grRules=grRules(1:counter); +0674 rxnMiriams=rxnMiriams(1:counter); +0675 reactionReversibility=reactionReversibility(1:counter); +0676 reactionUB=reactionUB(1:counter); +0677 reactionLB=reactionLB(1:counter); +0678 reactionObjective=reactionObjective(1:counter); +0679 S=S(:,1:counter); +0680 +0681 model.name=modelSBML.name; +0682 model.id=modelSBML.id; +0683 model.rxns=reactionIDs; +0684 model.mets=metaboliteIDs; +0685 model.S=sparse(S); +0686 model.lb=reactionLB; +0687 model.ub=reactionUB; +0688 model.rev=reactionReversibility; +0689 model.c=reactionObjective; +0690 model.b=zeros(numel(metaboliteIDs),1); +0691 model.comps=compartmentIDs; +0692 model.compNames=compartmentNames; +0693 model.rxnConfidenceScores=rxnconfidencescores; +0694 model.rxnReferences=rxnreferences; +0695 model.rxnNotes=rxnnotes; +0696 +0697 %Load annotation if available. If there are several authors, only the first +0698 %author credentials are imported +0699 if isfield(modelSBML,'annotation') +0700 endString='</'; +0701 I=strfind(modelSBML.annotation,endString); +0702 J=strfind(modelSBML.annotation,'<vCard:Family>'); +0703 if any(J) +0704 model.annotation.familyName=modelSBML.annotation(J(1)+14:I(find(I>J(1),1))-1); +0705 end +0706 J=strfind(modelSBML.annotation,'<vCard:Given>'); +0707 if any(J) +0708 model.annotation.givenName=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); +0709 end +0710 J=strfind(modelSBML.annotation,'<vCard:EMAIL>'); +0711 if any(J) +0712 model.annotation.email=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); +0713 end +0714 J=strfind(modelSBML.annotation,'<vCard:Orgname>'); +0715 if any(J) +0716 model.annotation.organization=modelSBML.annotation(J(1)+15:I(find(I>J(1),1))-1); +0717 end +0718 endString='"/>'; +0719 I=strfind(modelSBML.annotation,endString); +0720 if strfind(modelSBML.annotation,'"urn:miriam:') +0721 J=strfind(modelSBML.annotation,'"urn:miriam:'); +0722 if any(J) +0723 model.annotation.taxonomy=modelSBML.annotation(J+12:I(find(I>J,1))-1); +0724 end +0725 else +0726 J=strfind(modelSBML.annotation,'"http://identifiers.org/'); +0727 if any(J) +0728 I = I(find(I>J,1))-1; +0729 J = J(find(J<I,1))+24; +0730 model.annotation.taxonomy=modelSBML.annotation(J:I); +0731 else +0732 J=strfind(modelSBML.annotation,'"https://identifiers.org/'); +0733 if any(J) +0734 I = I(find(I>J,1))-1; +0735 J = J(find(J<I,1))+25; +0736 model.annotation.taxonomy=modelSBML.annotation(J:I); +0737 end +0738 end +0739 end +0740 end +0741 if isfield(modelSBML,'notes') +0742 startString=strfind(modelSBML.notes,'xhtml">'); +0743 endString=strfind(modelSBML.notes,'</body>'); +0744 if any(startString) && any(endString) +0745 model.annotation.note=modelSBML.notes(startString+7:endString-1); +0746 model.annotation.note=regexprep(model.annotation.note,'<p>|</p>',''); +0747 model.annotation.note=strtrim(model.annotation.note); +0748 if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') +0749 model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel +0750 end +0751 end +0752 end +0753 +0754 if any(~cellfun(@isempty,compartmentOutside)) +0755 model.compOutside=compartmentOutside; +0756 end +0757 +0758 model.rxnNames=reactionNames; +0759 model.metNames=metaboliteNames; +0760 +0761 %Match the compartments for metabolites +0762 [~, J]=ismember(metaboliteCompartments,model.comps); +0763 model.metComps=J; +0764 +0765 %If any genes have been loaded (only for the new format) +0766 if ~isempty(geneNames) +0767 %In some rare cases geneNames may not necessarily be used in grRules. +0768 %That is true for Yeast 7.6. It's therefore important to change gene +0769 %systematic names to geneIDs in sophisticated way. Gene systematic +0770 %names are not unique, since exactly the same name may be in different +0771 %compartments +0772 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) +0773 geneShortNames=geneNames; +0774 %geneShortNames contain compartments as well, so these are removed +0775 geneShortNames=regexprep(geneShortNames,' \[.+$',''); +0776 %grRules obtained from modifier fields contain geneNames. These are +0777 %changed into geneIDs. grRulesFromModifier is a good way to have +0778 %geneIDs and rxns association when it's important to resolve +0779 %systematic name ambiguities +0780 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0781 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0782 +0783 %Yeast 7.6 contains several metabolites, which were used in gene +0784 %associations. For that reason, the list of species ID is created +0785 %and we then check whether any of them have kegg.genes annotation +0786 %thereby obtaining systematic gene names +0787 geneShortNames=vertcat(geneShortNames,metaboliteNames); +0788 geneIDs=vertcat(geneIDs,metaboliteIDs); +0789 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); +0790 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); +0791 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); +0792 +0793 %Now we retain information for only these entries, which have +0794 %kegg.genes annotation +0795 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); +0796 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); +0797 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); +0798 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); +0799 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); +0800 %Now we reorder geneIDs and geneSystNames by geneSystNames string +0801 %length +0802 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later +0803 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); +0804 geneIDs = geneIDs(Indx); +0805 geneSystNames = geneSystNames(Indx); +0806 for i=1:numel(geneSystNames) +0807 for j=1:numel(grRules) +0808 if strfind(grRules{j},geneSystNames{i}) +0809 if ~isempty(grRules{j}) +0810 if sum(ismember(geneSystNames,geneSystNames{i}))==1 +0811 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); +0812 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 +0813 counter=0; +0814 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); +0815 for k=1:numel(ovrlpIDs) +0816 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) +0817 counter=counter+1; +0818 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); +0819 end +0820 if counter>1 +0821 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; +0822 dispEM(EM); +0823 end +0824 end +0825 end +0826 end +0827 end +0828 end +0829 end +0830 end +0831 model.genes=geneNames; +0832 model.grRules=grRules; +0833 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0834 model.grRules = grRules; +0835 model.rxnGeneMat = rxnGeneMat; +0836 +0837 %Match the compartments for genes +0838 [~, J]=ismember(geneCompartments,model.comps); +0839 model.geneComps=J; +0840 else +0841 if ~all(cellfun(@isempty,grRules)) +0842 %If fbc_geneProduct exists, follow the specified gene order, such +0843 %that matching geneShortNames in function below will work +0844 if isfield(modelSBML,'fbc_geneProduct') +0845 genes={modelSBML.fbc_geneProduct.fbc_id}; +0846 +0847 %Get gene Miriams if they were not retrieved above (this occurs +0848 %when genes are stored as fbc_geneProduct instead of species) +0849 if isempty(geneMiriams) +0850 geneMiriams = cell(numel(genes),1); +0851 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm])) +0852 %If all the SBO terms are identical, don't add them to geneMiriams +0853 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); +0854 end +0855 for i = 1:numel(genes) +0856 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); +0857 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1) +0858 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); +0859 end 0860 end -0861 for i = 1:numel(genes) -0862 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); -0863 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1) -0864 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); -0865 end -0866 end -0867 end -0868 proteins={modelSBML.fbc_geneProduct.fbc_name}; -0869 else -0870 genes=getGenesFromGrRules(grRules); -0871 end -0872 model.genes=genes; -0873 model.grRules=grRules; -0874 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0875 model.grRules = grRules; -0876 model.rxnGeneMat = rxnGeneMat; -0877 end -0878 end -0879 -0880 if all(cellfun(@isempty,geneShortNames)) -0881 if isfield(modelSBML,'fbc_geneProduct') -0882 for i=1:numel(genes) -0883 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) -0884 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; -0885 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) -0886 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; -0887 else -0888 geneShortNames{i,1}=''; -0889 end -0890 end -0891 end -0892 end -0893 -0894 %If any InChIs have been loaded -0895 if any(~cellfun(@isempty,metaboliteInChI)) -0896 model.inchis=metaboliteInChI; -0897 end -0898 -0899 %If any formulas have been loaded -0900 if any(~cellfun(@isempty,metaboliteFormula)) -0901 model.metFormulas=metaboliteFormula; -0902 end -0903 -0904 %If any charges have been loaded -0905 if ~isempty(metaboliteCharges) -0906 model.metCharges=metaboliteCharges; -0907 end -0908 -0909 %If any gene short names have been loaded -0910 if any(~cellfun(@isempty,geneShortNames)) -0911 model.geneShortNames=geneShortNames; -0912 end -0913 -0914 %If any Miriam strings for compartments have been loaded -0915 if any(~cellfun(@isempty,compartmentMiriams)) -0916 model.compMiriams=compartmentMiriams; -0917 end -0918 -0919 %If any Miriam strings for metabolites have been loaded -0920 if any(~cellfun(@isempty,metaboliteMiriams)) -0921 model.metMiriams=metaboliteMiriams; -0922 end -0923 -0924 %If any subsystems have been loaded -0925 if any(~cellfun(@isempty,subsystems)) -0926 model.subSystems=subsystems; -0927 end -0928 if any(rxnComps) -0929 if all(rxnComps) -0930 model.rxnComps=rxnComps; -0931 else -0932 if supressWarnings==false -0933 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; -0934 dispEM(EM,false,model.rxns(rxnComps==0)); -0935 end -0936 end -0937 end -0938 -0939 %If any ec-codes have been loaded -0940 if any(~cellfun(@isempty,eccodes)) -0941 model.eccodes=eccodes; -0942 end -0943 -0944 %If any Miriam strings for reactions have been loaded -0945 if any(~cellfun(@isempty,rxnMiriams)) -0946 model.rxnMiriams=rxnMiriams; -0947 end -0948 -0949 %If any Miriam strings for genes have been loaded -0950 if any(~cellfun(@isempty,geneMiriams)) -0951 model.geneMiriams=geneMiriams; +0861 end +0862 proteins={modelSBML.fbc_geneProduct.fbc_name}; +0863 else +0864 genes=getGenesFromGrRules(grRules); +0865 end +0866 model.genes=genes; +0867 model.grRules=grRules; +0868 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0869 model.grRules = grRules; +0870 model.rxnGeneMat = rxnGeneMat; +0871 end +0872 end +0873 +0874 if all(cellfun(@isempty,geneShortNames)) +0875 if isfield(modelSBML,'fbc_geneProduct') +0876 for i=1:numel(genes) +0877 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) +0878 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; +0879 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) +0880 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; +0881 else +0882 geneShortNames{i,1}=''; +0883 end +0884 end +0885 end +0886 end +0887 +0888 %If any InChIs have been loaded +0889 if any(~cellfun(@isempty,metaboliteInChI)) +0890 model.inchis=metaboliteInChI; +0891 end +0892 +0893 %If any formulas have been loaded +0894 if any(~cellfun(@isempty,metaboliteFormula)) +0895 model.metFormulas=metaboliteFormula; +0896 end +0897 +0898 %If any charges have been loaded +0899 if ~isempty(metaboliteCharges) +0900 model.metCharges=metaboliteCharges; +0901 end +0902 +0903 %If any gene short names have been loaded +0904 if any(~cellfun(@isempty,geneShortNames)) +0905 model.geneShortNames=geneShortNames; +0906 end +0907 +0908 %If any Miriam strings for compartments have been loaded +0909 if any(~cellfun(@isempty,compartmentMiriams)) +0910 model.compMiriams=compartmentMiriams; +0911 end +0912 +0913 %If any Miriam strings for metabolites have been loaded +0914 if any(~cellfun(@isempty,metaboliteMiriams)) +0915 model.metMiriams=metaboliteMiriams; +0916 end +0917 +0918 %If any subsystems have been loaded +0919 if any(~cellfun(@isempty,subsystems)) +0920 model.subSystems=subsystems; +0921 end +0922 if any(rxnComps) +0923 if all(rxnComps) +0924 model.rxnComps=rxnComps; +0925 else +0926 if supressWarnings==false +0927 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; +0928 dispEM(EM,false,model.rxns(rxnComps==0)); +0929 end +0930 end +0931 end +0932 +0933 %If any ec-codes have been loaded +0934 if any(~cellfun(@isempty,eccodes)) +0935 model.eccodes=eccodes; +0936 end +0937 +0938 %If any Miriam strings for reactions have been loaded +0939 if any(~cellfun(@isempty,rxnMiriams)) +0940 model.rxnMiriams=rxnMiriams; +0941 end +0942 +0943 %If any Miriam strings for genes have been loaded +0944 if any(~cellfun(@isempty,geneMiriams)) +0945 model.geneMiriams=geneMiriams; +0946 end +0947 +0948 %If any protein strings have been loaded +0949 if any(~cellfun(@isempty,proteins)) +0950 proteins = reshape(proteins,[],1); +0951 model.proteins=proteins; 0952 end 0953 -0954 %If any protein strings have been loaded -0955 if any(~cellfun(@isempty,proteins)) -0956 proteins = reshape(proteins,[],1); -0957 model.proteins=proteins; -0958 end -0959 -0960 model.unconstrained=metaboliteUnconstrained; -0961 -0962 %Convert SBML IDs back into their original strings. Here we are using part -0963 %from convertSBMLID, originating from the COBRA Toolbox -0964 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); -0965 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); -0966 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); -0967 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); -0968 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); -0969 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); -0970 -0971 if removePrefix -0972 [model, hasChanged]=removeIdentifierPrefix(model); -0973 dispEM(['The following fields have prefixes removed from all entries. '... -0974 'If this is undesired, run importModel with removePrefix as false. Example: '... -0975 'importModel(''filename.xml'',[],false);'],false,hasChanged) -0976 end -0977 -0978 %Remove unused fields -0979 if isempty(model.annotation) -0980 model=rmfield(model,'annotation'); +0954 model.unconstrained=metaboliteUnconstrained; +0955 +0956 %Convert SBML IDs back into their original strings. Here we are using part +0957 %from convertSBMLID, originating from the COBRA Toolbox +0958 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); +0959 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); +0960 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); +0961 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); +0962 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); +0963 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); +0964 +0965 if removePrefix +0966 [model, hasChanged]=removeIdentifierPrefix(model); +0967 dispEM(['The following fields have prefixes removed from all entries. '... +0968 'If this is undesired, run importModel with removePrefix as false. Example: '... +0969 'importModel(''filename.xml'',[],false);'],false,hasChanged) +0970 end +0971 +0972 %Remove unused fields +0973 if isempty(model.annotation) +0974 model=rmfield(model,'annotation'); +0975 end +0976 if isempty(model.compOutside) +0977 model=rmfield(model,'compOutside'); +0978 end +0979 if isempty(model.compMiriams) +0980 model=rmfield(model,'compMiriams'); 0981 end -0982 if isempty(model.compOutside) -0983 model=rmfield(model,'compOutside'); +0982 if isempty(model.rxnComps) +0983 model=rmfield(model,'rxnComps'); 0984 end -0985 if isempty(model.compMiriams) -0986 model=rmfield(model,'compMiriams'); +0985 if isempty(model.grRules) +0986 model=rmfield(model,'grRules'); 0987 end -0988 if isempty(model.rxnComps) -0989 model=rmfield(model,'rxnComps'); +0988 if isempty(model.rxnGeneMat) +0989 model=rmfield(model,'rxnGeneMat'); 0990 end -0991 if isempty(model.grRules) -0992 model=rmfield(model,'grRules'); -0993 end -0994 if isempty(model.rxnGeneMat) -0995 model=rmfield(model,'rxnGeneMat'); -0996 end -0997 if isempty(model.subSystems) -0998 model=rmfield(model,'subSystems'); -0999 else -1000 model.subSystems(cellfun(@isempty,subsystems))={{''}}; +0991 if isempty(model.subSystems) +0992 model=rmfield(model,'subSystems'); +0993 else +0994 model.subSystems(cellfun(@isempty,subsystems))={{''}}; +0995 end +0996 if isempty(model.eccodes) +0997 model=rmfield(model,'eccodes'); +0998 end +0999 if isempty(model.rxnMiriams) +1000 model=rmfield(model,'rxnMiriams'); 1001 end -1002 if isempty(model.eccodes) -1003 model=rmfield(model,'eccodes'); +1002 if cellfun(@isempty,model.rxnNotes) +1003 model=rmfield(model,'rxnNotes'); 1004 end -1005 if isempty(model.rxnMiriams) -1006 model=rmfield(model,'rxnMiriams'); +1005 if cellfun(@isempty,model.rxnReferences) +1006 model=rmfield(model,'rxnReferences'); 1007 end -1008 if cellfun(@isempty,model.rxnNotes) -1009 model=rmfield(model,'rxnNotes'); +1008 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) +1009 model=rmfield(model,'rxnConfidenceScores'); 1010 end -1011 if cellfun(@isempty,model.rxnReferences) -1012 model=rmfield(model,'rxnReferences'); -1013 end -1014 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) -1015 model=rmfield(model,'rxnConfidenceScores'); -1016 end -1017 if isempty(model.genes) -1018 model=rmfield(model,'genes'); -1019 elseif isrow(model.genes) -1020 model.genes=transpose(model.genes); +1011 if isempty(model.genes) +1012 model=rmfield(model,'genes'); +1013 elseif isrow(model.genes) +1014 model.genes=transpose(model.genes); +1015 end +1016 if isempty(model.geneComps) +1017 model=rmfield(model,'geneComps'); +1018 end +1019 if isempty(model.geneMiriams) +1020 model=rmfield(model,'geneMiriams'); 1021 end -1022 if isempty(model.geneComps) -1023 model=rmfield(model,'geneComps'); +1022 if isempty(model.geneShortNames) +1023 model=rmfield(model,'geneShortNames'); 1024 end -1025 if isempty(model.geneMiriams) -1026 model=rmfield(model,'geneMiriams'); +1025 if isempty(model.proteins) +1026 model=rmfield(model,'proteins'); 1027 end -1028 if isempty(model.geneShortNames) -1029 model=rmfield(model,'geneShortNames'); +1028 if isempty(model.inchis) +1029 model=rmfield(model,'inchis'); 1030 end -1031 if isempty(model.proteins) -1032 model=rmfield(model,'proteins'); +1031 if isempty(model.metFormulas) +1032 model=rmfield(model,'metFormulas'); 1033 end -1034 if isempty(model.inchis) -1035 model=rmfield(model,'inchis'); +1034 if isempty(model.metMiriams) +1035 model=rmfield(model,'metMiriams'); 1036 end -1037 if isempty(model.metFormulas) -1038 model=rmfield(model,'metFormulas'); +1037 if ~any(model.metCharges) +1038 model=rmfield(model,'metCharges'); 1039 end -1040 if isempty(model.metMiriams) -1041 model=rmfield(model,'metMiriams'); -1042 end -1043 if ~any(model.metCharges) -1044 model=rmfield(model,'metCharges'); -1045 end -1046 -1047 %This just removes the grRules if no genes have been loaded -1048 if ~isfield(model,'genes') && isfield(model,'grRules') -1049 model=rmfield(model,'grRules'); -1050 end -1051 -1052 %Print warnings about bad structure -1053 if supressWarnings==false -1054 checkModelStruct(model,false); -1055 end -1056 -1057 if removeExcMets==true -1058 model=simplifyModel(model); -1059 end -1060 end +1040 +1041 %This just removes the grRules if no genes have been loaded +1042 if ~isfield(model,'genes') && isfield(model,'grRules') +1043 model=rmfield(model,'grRules'); +1044 end +1045 +1046 %Print warnings about bad structure +1047 if supressWarnings==false +1048 checkModelStruct(model,false); +1049 end +1050 +1051 if removeExcMets==true +1052 model=simplifyModel(model); +1053 end +1054 end +1055 +1056 function fieldContent=parseNote(searchString,fieldName) +1057 %The function obtains the particular information from 'notes' field, using +1058 %fieldName as the dummy string +1059 +1060 fieldContent=''; 1061 -1062 function fieldContent=parseNote(searchString,fieldName) -1063 %The function obtains the particular information from 'notes' field, using -1064 %fieldName as the dummy string -1065 -1066 fieldContent=''; -1067 -1068 if strfind(searchString,fieldName) -1069 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); -1070 targetString=regexprep(targetString,'<p>|</p>',''); -1071 targetString=regexprep(targetString,[fieldName, ':'],''); -1072 for i=1:numel(targetString) -1073 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1074 end -1075 fieldContent=regexprep(fieldContent,'^;|;$',''); -1076 else -1077 fieldContent=''; -1078 end -1079 end -1080 -1081 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) -1082 -1083 fieldContent=''; -1084 -1085 %Removing whitespace characters from the ending strings, which may occur in -1086 %several cases -1087 searchString=regexprep(searchString,'" />','"/>'); -1088 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); -1089 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1090 targetString=regexprep(targetString,startString,''); -1091 targetString=regexprep(targetString,[fieldName midString],''); -1092 -1093 for i=1:numel(targetString) -1094 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1095 end +1062 if strfind(searchString,fieldName) +1063 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); +1064 targetString=regexprep(targetString,'<p>|</p>',''); +1065 targetString=regexprep(targetString,[fieldName, ':'],''); +1066 for i=1:numel(targetString) +1067 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1068 end +1069 fieldContent=regexprep(fieldContent,'^;|;$',''); +1070 else +1071 fieldContent=''; +1072 end +1073 end +1074 +1075 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) +1076 +1077 fieldContent=''; +1078 +1079 %Removing whitespace characters from the ending strings, which may occur in +1080 %several cases +1081 searchString=regexprep(searchString,'" />','"/>'); +1082 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); +1083 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1084 targetString=regexprep(targetString,startString,''); +1085 targetString=regexprep(targetString,[fieldName midString],''); +1086 +1087 for i=1:numel(targetString) +1088 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1089 end +1090 +1091 fieldContent=regexprep(fieldContent,'^;|;$',''); +1092 end +1093 +1094 function miriamStruct=parseMiriam(searchString) +1095 %Generates miriam structure from annotation field 1096 -1097 fieldContent=regexprep(fieldContent,'^;|;$',''); -1098 end -1099 -1100 function miriamStruct=parseMiriam(searchString) -1101 %Generates miriam structure from annotation field -1102 -1103 %Finding whether miriams are written in the old or the new way -1104 if strfind(searchString,'urn:miriam:') -1105 startString='urn:miriam:'; -1106 elseif strfind(searchString,'http://identifiers.org/') -1107 startString='http://identifiers.org/'; -1108 elseif strfind(searchString,'https://identifiers.org/') -1109 startString='https://identifiers.org/'; -1110 else -1111 miriamStruct=[]; -1112 return; -1113 end -1114 -1115 miriamStruct=[]; -1116 -1117 searchString=regexprep(searchString,'" />','"/>'); -1118 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); -1119 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1120 targetString=regexprep(targetString,startString,''); -1121 -1122 fwdslash = contains(targetString,'/'); -1123 midString = cell(numel(targetString),1); -1124 midString(fwdslash) = {'/'}; -1125 midString(~fwdslash) = {':'}; -1126 -1127 counter=0; -1128 for i=1:numel(targetString) -1129 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) -1130 counter=counter+1; -1131 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); -1132 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); -1133 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); -1134 end -1135 end -1136 end -1137 -1138 function miriam = addSBOtoMiriam(miriam,sboTerm) -1139 %Appends SBO term to miriam structure -1140 -1141 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format -1142 if isempty(miriam) -1143 miriam.name = {'sbo'}; -1144 miriam.value = sboTerm; -1145 elseif any(strcmp('sbo',miriam.name)) -1146 currSbo = strcmp('sbo',miriam.name); -1147 miriam.value(currSbo) = sboTerm; -1148 else -1149 miriam.name(end+1) = {'sbo'}; -1150 miriam.value(end+1) = sboTerm; -1151 end -1152 end +1097 %Finding whether miriams are written in the old or the new way +1098 if strfind(searchString,'urn:miriam:') +1099 startString='urn:miriam:'; +1100 elseif strfind(searchString,'http://identifiers.org/') +1101 startString='http://identifiers.org/'; +1102 elseif strfind(searchString,'https://identifiers.org/') +1103 startString='https://identifiers.org/'; +1104 else +1105 miriamStruct=[]; +1106 return; +1107 end +1108 +1109 miriamStruct=[]; +1110 +1111 searchString=regexprep(searchString,'" />','"/>'); +1112 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); +1113 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1114 targetString=regexprep(targetString,startString,''); +1115 +1116 fwdslash = contains(targetString,'/'); +1117 midString = cell(numel(targetString),1); +1118 midString(fwdslash) = {'/'}; +1119 midString(~fwdslash) = {':'}; +1120 +1121 counter=0; +1122 for i=1:numel(targetString) +1123 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) +1124 counter=counter+1; +1125 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); +1126 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); +1127 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); +1128 end +1129 end +1130 end +1131 +1132 function miriam = addSBOtoMiriam(miriam,sboTerm) +1133 %Appends SBO term to miriam structure +1134 +1135 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format +1136 if isempty(miriam) +1137 miriam.name = {'sbo'}; +1138 miriam.value = sboTerm; +1139 elseif any(strcmp('sbo',miriam.name)) +1140 currSbo = strcmp('sbo',miriam.name); +1141 miriam.value(currSbo) = sboTerm; +1142 else +1143 miriam.name(end+1) = {'sbo'}; +1144 miriam.value(end+1) = sboTerm; +1145 end +1146 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/solver/solveLP.html b/doc/solver/solveLP.html index 4532a550..8c462e72 100644 --- a/doc/solver/solveLP.html +++ b/doc/solver/solveLP.html @@ -43,7 +43,7 @@

    DESCRIPTION ^SOURCE CODE ^% interpret. Note that this optimization can be very 0016 % slow 0017 % (optional, default 0) -0018 % params parameter structure as used by getMILPParams (optional) +0018 % params *obsolete option* 0019 % hsSol hot-start solution for the LP solver. This can 0020 % significantly speed up the process if many similar 0021 % optimization problems are solved iteratively. Only used diff --git a/io/importModel.m b/io/importModel.m index e8a412d4..a441d94f 100755 --- a/io/importModel.m +++ b/io/importModel.m @@ -155,7 +155,7 @@ compartmentOutside=cell(numel(modelSBML.compartment),1); compartmentMiriams=cell(numel(modelSBML.compartment),1); -if isfield(modelSBML.compartment,'sboTerm') && numel(unique([modelSBML.compartment.sboTerm])) == 1 +if isfield(modelSBML.compartment,'sboTerm') && isscalar(unique([modelSBML.compartment.sboTerm])) %If all the SBO terms are identical, don't add them to compMiriams modelSBML.compartment = rmfield(modelSBML.compartment,'sboTerm'); end @@ -291,7 +291,7 @@ metaboliteFormula{numel(metaboliteFormula)+1,1}=... formula(compositionIndexes(1)+1:compositionIndexes(2)-1); else - if numel(compositionIndexes)==1 + if isscalar(compositionIndexes) %Probably a simple molecule which can have only %one conformation metaboliteFormula{numel(metaboliteFormula)+1,1}=... @@ -422,7 +422,7 @@ parameter.value={modelSBML.parameter(:).value}'; end -if isfield(modelSBML.reaction,'sboTerm') && numel(unique([modelSBML.reaction.sboTerm])) == 1 +if isfield(modelSBML.reaction,'sboTerm') && isscalar(unique([modelSBML.reaction.sboTerm])) %If all the SBO terms are identical, don't add them to rxnMiriams modelSBML.reaction = rmfield(modelSBML.reaction,'sboTerm'); end @@ -433,7 +433,7 @@ %so, then jump to the next reaction. This is because I get the genes %for complexes from the names and not from the reactions that create %them. This only applies to the non-COBRA format - if numel(modelSBML.reaction(i).product)==1 + if isscalar(modelSBML.reaction(i).product) if length(modelSBML.reaction(i).product(1).species)>=3 if strcmp(modelSBML.reaction(i).product(1).species(1:3),'Cx_')==true continue; @@ -455,19 +455,9 @@ if isfield(modelSBML.reaction(i),'fbc_lowerFluxBound') lb=modelSBML.reaction(i).fbc_lowerFluxBound; ub=modelSBML.reaction(i).fbc_upperFluxBound; - for n=1:numel(parameter.value) - lb=regexprep(lb,parameter.name(n),num2str(parameter.value{n})); - ub=regexprep(ub,parameter.name(n),num2str(parameter.value{n})); - end - if isempty(lb) - lb='-Inf'; - end - if isempty(ub) - ub='Inf'; - end - reactionLB(counter)=str2num(lb); - reactionUB(counter)=str2num(ub); - %The order of these parameters should not be hard coded + [~,fluxBoundIdx] = ismember({lb,ub},parameter.name); + reactionLB(counter) = parameter.value{fluxBoundIdx(1)}; + reactionUB(counter) = parameter.value{fluxBoundIdx(2)}; elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter') reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value; reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value; @@ -541,7 +531,7 @@ %genes from the name of the complex (not the %reaction that creates it) index=find(strcmp(modifier,complexIDs)); - if numel(index)==1 + if isscalar(index) if ~isempty(rules) rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')']; else @@ -648,7 +638,7 @@ if strcmp(obj,modelSBML.fbc_objective(i).fbc_id) if ~isempty(modelSBML.fbc_objective(i).fbc_fluxObjective) rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction; - idx=find(ismember(reactionIDs,rxn)); + idx=ismember(reactionIDs,rxn); reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient; end end @@ -735,11 +725,15 @@ else J=strfind(modelSBML.annotation,'"http://identifiers.org/'); if any(J) - model.annotation.taxonomy=modelSBML.annotation(J+24:I(find(I>J,1))-1); + I = I(find(I>J,1))-1; + J = J(find(JJ,1))-1); + I = I(find(I>J,1))-1; + J = J(find(J