removeBadRxns Iteratively removes reactions which enable production/consumption of some metabolite without any uptake/excretion model a model structure. For the intented function, the model shouldn't allow for any uptake/excretion. The easiest way to achieve this is to import the model using importModel('filename',false) rxnRules 1: only remove reactions which are unbalanced 2: also remove reactions which couldn't be checked for mass balancing 3: all reactions can be removed (opt, default 1) ignoreMets either a cell array of metabolite IDs, a logical vector with the same number of elements as metabolites in the model, of a vector of indexes for metabolites to exclude from this analysis (opt, default []) isNames true if the supplied mets represent metabolite names (as opposed to IDs). This is a way to delete metabolites in several compartments at once without knowing the exact IDs. This only works if ignoreMets is a cell array (opt, default false) balanceElements a cell array with the elements for which to balance the reactions. May contain any combination of 'C','P','S','N','R','O','H' (where R stands for R group) (opt, default {'C';'P';'S';'N';'O'} refModel a reference model which can be used to ensure that the resulting model is still functional. The intended use is that the reference model is a copy of model, but with uptake/excretion allowed and some objectives (such as production of biomass) constrained to a non-zero flux. Before a reaction is removed from "model" the function first checks that the same deletion in "refModel" doesn't render the problem unfeasible (opt) 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 not possible using the default constraints (opt, default false) printReport true if a report should be printed (opt, default false) newModel a model structure after the problematic reactions have been deleted removedRxns a cell array with the reactions that were removed The purpose of this function is to remove reactions which enable production/consumption of metabolites even when exchange reactions aren't used. Many models, especially if they are automatically inferred from databases, will have unbalanced reactions which allow for net-production/consumption of metabolites without any consumption/excretion. A common reason for this is when general compounds have different meaning in different reactions (as DNA has in these two reactions). dATP + dGTP + dCTP + dTTP <=> DNA + 4 PPi 0.25 dATP + 0.25 dGTP + 0.25 dCTP + 0.25 dTTP <=> DNA + PPi Reactions that are problematic like this are always elementally unbalanced, but it is not always that you would like to exclude all unbalanced reactions from your model. This function tries to remove as few problematic reactions as possible so that the model cannot produce/consume anything from nothing. This is done by repeatedly calling makeSomething/consumeSomething, checking if any of the involved reactions are elementally unbalanced, remove one of them, and then iterating until no metabolites can be produced/consumed. makeSomething is called before consumeSomething. Usage: [newModel removedRxns]=removeBadRxns(model,rxnRules,... ignoreMets,isNames,refModel,ignoreIntBounds,printReport) Rasmus Agren, 2013-05-16
0001 function [newModel removedRxns]=removeBadRxns(model,rxnRules,ignoreMets,isNames,balanceElements,refModel,ignoreIntBounds,printReport) 0002 % removeBadRxns 0003 % Iteratively removes reactions which enable production/consumption of some 0004 % metabolite without any uptake/excretion 0005 % 0006 % model a model structure. For the intented function, 0007 % the model shouldn't allow for any uptake/excretion. 0008 % The easiest way to achieve this is to import the 0009 % model using importModel('filename',false) 0010 % rxnRules 1: only remove reactions which are unbalanced 0011 % 2: also remove reactions which couldn't be checked for 0012 % mass balancing 0013 % 3: all reactions can be removed 0014 % (opt, default 1) 0015 % ignoreMets either a cell array of metabolite IDs, a logical vector 0016 % with the same number of elements as metabolites in the model, 0017 % of a vector of indexes for metabolites to exclude from 0018 % this analysis (opt, default []) 0019 % isNames true if the supplied mets represent metabolite names 0020 % (as opposed to IDs). This is a way to delete 0021 % metabolites in several compartments at once without 0022 % knowing the exact IDs. This only works if ignoreMets 0023 % is a cell array (opt, default false) 0024 % balanceElements a cell array with the elements for which to 0025 % balance the reactions. May contain any 0026 % combination of 'C','P','S','N','R','O','H' 0027 % (where R stands for R group) (opt, default 0028 % {'C';'P';'S';'N';'O'} 0029 % refModel a reference model which can be used to ensure 0030 % that the resulting model is still functional. 0031 % The intended use is that the reference model is 0032 % a copy of model, but with uptake/excretion allowed and 0033 % some objectives (such as production of biomass) 0034 % constrained to a non-zero flux. Before a 0035 % reaction is removed from "model" the function first 0036 % checks that the same deletion in "refModel" 0037 % doesn't render the problem unfeasible (opt) 0038 % ignoreIntBounds true if internal bounds (including reversibility) 0039 % should be ignored. Exchange reactions are not affected. 0040 % This can be used to find unbalanced solutions which are 0041 % not possible using the default constraints (opt, 0042 % default false) 0043 % printReport true if a report should be printed (opt, 0044 % default false) 0045 % 0046 % newModel a model structure after the problematic 0047 % reactions have been deleted 0048 % removedRxns a cell array with the reactions that were 0049 % removed 0050 % 0051 % The purpose of this function is to remove reactions which enable 0052 % production/consumption of metabolites even when exchange reactions aren't used. 0053 % Many models, especially if they are automatically inferred from 0054 % databases, will have unbalanced reactions which allow for 0055 % net-production/consumption of metabolites without any consumption/excretion. 0056 % A common reason for this is when general compounds have different meaning 0057 % in different reactions (as DNA has in these two reactions). 0058 % dATP + dGTP + dCTP + dTTP <=> DNA + 4 PPi 0059 % 0.25 dATP + 0.25 dGTP + 0.25 dCTP + 0.25 dTTP <=> DNA + PPi 0060 % Reactions that are problematic like this are always elementally 0061 % unbalanced, but it is not always that you would like to exclude all 0062 % unbalanced reactions from your model. 0063 % This function tries to remove as few problematic reactions as possible 0064 % so that the model cannot produce/consume anything from nothing. This is done by 0065 % repeatedly calling makeSomething/consumeSomething, checking if any of 0066 % the involved reactions are elementally unbalanced, remove one of them, 0067 % and then iterating until no metabolites can be produced/consumed. 0068 % makeSomething is called before consumeSomething. 0069 % 0070 % Usage: [newModel removedRxns]=removeBadRxns(model,rxnRules,... 0071 % ignoreMets,isNames,refModel,ignoreIntBounds,printReport) 0072 % 0073 % Rasmus Agren, 2013-05-16 0074 % 0075 0076 if nargin<2 0077 rxnRules=1; 0078 end 0079 if nargin<3 0080 ignoreMets=[]; 0081 end 0082 if nargin<4 0083 isNames=false; 0084 end 0085 if nargin<5 0086 balanceElements={'C';'P';'S';'N';'O'}; 0087 end 0088 if nargin<6 0089 refModel=[]; 0090 else 0091 if ~isempty(refModel) 0092 fprintf('WARNING: The feature to supply a reference model is currently not supported\n'); 0093 end 0094 end 0095 if nargin<7 0096 ignoreIntBounds=false; 0097 end 0098 if nargin<8 0099 printReport=false; 0100 end 0101 0102 %Check if the model has open exchange reactions and print a warning in that case 0103 if ~isfield(model,'unconstrained') 0104 [crap I]=getExchangeRxns(model); 0105 if any(I) 0106 if any(model.lb(I)~=0) || any(model.ub(I)~=0) 0107 fprintf('WARNING: The model contains open exchange reactions. This is not the intended use of this function. Consider importing your model using importModel(filename,false).\n'); 0108 end 0109 end 0110 end 0111 0112 %Check that the model is feasible 0113 sol=solveLP(model); 0114 if isempty(sol.f) 0115 throw(MException('','The model is not feasible. Consider removing lower bounds (such as ATP maintenance)')); 0116 end 0117 0118 %Check that the reference model is feasible 0119 if any(refModel) 0120 sol=solveLP(refModel); 0121 if isempty(sol.f) 0122 throw(MException('','The reference model is not feasible.')); 0123 end 0124 end 0125 0126 %Initialize some stuff 0127 removedRxns={}; 0128 balanceStructure=getElementalBalance(model); 0129 0130 %Check which elements to balance for 0131 if ~isempty(setdiff(upper(balanceElements),{'C','P','S','N','R','O','H'})) 0132 throw(MException('','Could not recognize all elements to balance for')); 0133 end 0134 bal=ismember({'C','P','S','N','R','O','H'},upper(balanceElements)); 0135 0136 %Main loop. First run for makeSomething, second for consumeSomething 0137 for i=1:2 0138 while 1 0139 %Make some metabolite using as few reactions as possible 0140 if i==1 0141 [solution metabolite]=makeSomething(model,ignoreMets,isNames,false,true,[],ignoreIntBounds); 0142 if ~isempty(solution) 0143 if printReport 0144 fprintf(['Can make: ' model.metNames{metabolite(1)} '\n']); 0145 end 0146 else 0147 %If no solution could be found, then finish 0148 break; 0149 end 0150 else 0151 [solution metabolite]=consumeSomething(model,ignoreMets,isNames,false,[],ignoreIntBounds); 0152 if ~isempty(solution) 0153 if printReport 0154 fprintf(['Can consume: ' model.metNames{metabolite(1)} '\n']); 0155 end 0156 else 0157 %If no solution could be found, then finish 0158 break; 0159 end 0160 end 0161 0162 %Find all reactions that are unbalanced and still carry flux 0163 I=find(abs(solution)>10^-8 & balanceStructure.balanceStatus>=0 & ~all(balanceStructure.leftComp(:,bal)==balanceStructure.rightComp(:,bal),2)); 0164 0165 %If there are unbalanced rxns then delete one of them and iterate 0166 if any(I) 0167 rxnToRemove=I(randsample(numel(I),1)); 0168 else 0169 %If there are no unbalanced rxns in the solution 0170 if rxnRules==1 0171 if i==1 0172 fprintf(['WARNING: No unbalanced reactions were found in the solution, but the model can still make "' model.metNames{metabolite} '". Aborting search. Consider setting rxnRules to 2 or 3 for a more exhaustive search\n']); 0173 else 0174 fprintf(['WARNING: No unbalanced reactions were found in the solution, but the model can still consume "' model.metNames{metabolite} '". Aborting search. Consider setting rxnRules to 2 or 3 for a more exhaustive search\n']); 0175 end 0176 break; 0177 else 0178 %Find reactions which are not checked for mass balancing, but 0179 %that still carry flux 0180 I=find(abs(solution)>10^-8 & balanceStructure.balanceStatus<0); 0181 0182 %If there are any such reactions, remove one of them and 0183 %iterate 0184 if any(I) 0185 rxnToRemove=I(randsample(numel(I),1)); 0186 else 0187 if rxnRules==2 0188 %This happens when all reactions used are balanced 0189 %according to the metabolite formulas. This cannot be, and 0190 %indicates that one or more of the formulas are wrong. 0191 %Print a warning and delete any reaction with flux 0192 if i==1 0193 fprintf(['WARNING: No unbalanced or unparsable reactions were found in the solution, but the model can still make "' model.metNames{metabolite} '". Aborting search. Consider setting rxnRules to 3 for a more exhaustive search\n']); 0194 else 0195 fprintf(['WARNING: No unbalanced or unparsable reactions were found in the solution, but the model can still consume "' model.metNames{metabolite} '". Aborting search. Consider setting rxnRules to 3 for a more exhaustive search\n']); 0196 end 0197 break; 0198 else 0199 if i==1 0200 fprintf(['WARNING: No unbalanced or unparsable reactions were found in the solution, but the model can still make "' model.metNames{metabolite} '". This indicates some error in the metabolite formulas. Removing a random reaction in the solution\n']); 0201 else 0202 fprintf(['WARNING: No unbalanced or unparsable reactions were found in the solution, but the model can still consume "' model.metNames{metabolite} '". This indicates some error in the metabolite formulas. Removing a random reaction in the solution\n']); 0203 end 0204 I=find(abs(solution)>10^-8); 0205 rxnToRemove=I(randsample(numel(I),1)); 0206 end 0207 end 0208 end 0209 end 0210 removedRxns=[removedRxns;model.rxns(rxnToRemove)]; 0211 if printReport 0212 fprintf(['\tRemoved: ' model.rxns{rxnToRemove} '\n']); 0213 end 0214 model=removeRxns(model,rxnToRemove); 0215 balanceStructure.balanceStatus(rxnToRemove)=[]; 0216 balanceStructure.leftComp(rxnToRemove,:)=[]; 0217 balanceStructure.rightComp(rxnToRemove,:)=[]; 0218 end 0219 end 0220 0221 newModel=model; 0222 end