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 the elements defined in parseFormulas (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-08-01
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 the elements defined in parseFormulas 0027 % (opt, default {'C';'P';'S';'N';'O'}) 0028 % refModel a reference model which can be used to ensure 0029 % that the resulting model is still functional. 0030 % The intended use is that the reference model is 0031 % a copy of model, but with uptake/excretion allowed and 0032 % some objectives (such as production of biomass) 0033 % constrained to a non-zero flux. Before a 0034 % reaction is removed from "model" the function first 0035 % checks that the same deletion in "refModel" 0036 % doesn't render the problem unfeasible (opt) 0037 % ignoreIntBounds true if internal bounds (including reversibility) 0038 % should be ignored. Exchange reactions are not affected. 0039 % This can be used to find unbalanced solutions which are 0040 % not possible using the default constraints (opt, 0041 % default false) 0042 % printReport true if a report should be printed (opt, 0043 % default false) 0044 % 0045 % newModel a model structure after the problematic 0046 % reactions have been deleted 0047 % removedRxns a cell array with the reactions that were 0048 % removed 0049 % 0050 % The purpose of this function is to remove reactions which enable 0051 % production/consumption of metabolites even when exchange reactions aren't used. 0052 % Many models, especially if they are automatically inferred from 0053 % databases, will have unbalanced reactions which allow for 0054 % net-production/consumption of metabolites without any consumption/excretion. 0055 % A common reason for this is when general compounds have different meaning 0056 % in different reactions (as DNA has in these two reactions). 0057 % dATP + dGTP + dCTP + dTTP <=> DNA + 4 PPi 0058 % 0.25 dATP + 0.25 dGTP + 0.25 dCTP + 0.25 dTTP <=> DNA + PPi 0059 % Reactions that are problematic like this are always elementally 0060 % unbalanced, but it is not always that you would like to exclude all 0061 % unbalanced reactions from your model. 0062 % This function tries to remove as few problematic reactions as possible 0063 % so that the model cannot produce/consume anything from nothing. This is done by 0064 % repeatedly calling makeSomething/consumeSomething, checking if any of 0065 % the involved reactions are elementally unbalanced, remove one of them, 0066 % and then iterating until no metabolites can be produced/consumed. 0067 % makeSomething is called before consumeSomething. 0068 % 0069 % Usage: [newModel removedRxns]=removeBadRxns(model,rxnRules,... 0070 % ignoreMets,isNames,refModel,ignoreIntBounds,printReport) 0071 % 0072 % Rasmus Agren, 2013-08-01 0073 % 0074 0075 if nargin<2 0076 rxnRules=1; 0077 end 0078 if nargin<3 0079 ignoreMets=[]; 0080 end 0081 if nargin<4 0082 isNames=false; 0083 end 0084 if nargin<5 0085 balanceElements={'C';'P';'S';'N';'O'}; 0086 end 0087 if nargin<6 0088 refModel=[]; 0089 else 0090 if ~isempty(refModel) 0091 dispEM('The feature to supply a reference model is currently not supported',false); 0092 end 0093 end 0094 if nargin<7 0095 ignoreIntBounds=false; 0096 end 0097 if nargin<8 0098 printReport=false; 0099 end 0100 0101 %Check if the model has open exchange reactions and print a warning in that case 0102 if ~isfield(model,'unconstrained') 0103 [crap I]=getExchangeRxns(model); 0104 if any(I) 0105 if any(model.lb(I)~=0) || any(model.ub(I)~=0) 0106 dispEM('The model contains open exchange reactions. This is not the intended use of this function. Consider importing your model using importModel(filename,false)',false); 0107 end 0108 end 0109 end 0110 0111 %Check that the model is feasible 0112 sol=solveLP(model); 0113 if isempty(sol.f) 0114 dispEM('The model is not feasible. Consider removing lower bounds (such as ATP maintenance)'); 0115 end 0116 0117 %Check that the reference model is feasible 0118 if any(refModel) 0119 sol=solveLP(refModel); 0120 if isempty(sol.f) 0121 dispEM('The reference model is not feasible'); 0122 end 0123 end 0124 0125 %Initialize some stuff 0126 removedRxns={}; 0127 balanceStructure=getElementalBalance(model); 0128 0129 %Check which elements to balance for 0130 if ~isempty(setdiff(balanceElements,balanceStructure.elements.abbrevs)) 0131 dispEM('Could not recognize all elements to balance for'); 0132 end 0133 bal=ismember(balanceStructure.elements.abbrevs,balanceElements); 0134 0135 %Main loop. First run for makeSomething, second for consumeSomething 0136 warned=false(2,1); %This is to prevent the same warning being printed multiple times if rxnRules==3 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 dispEM(['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'],false); 0173 else 0174 dispEM(['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'],false); 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 dispEM(['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'],false); 0194 else 0195 dispEM(['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'],false); 0196 end 0197 break; 0198 else 0199 if i==1 0200 if warned(1)==false 0201 dispEM(['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 random reactions in the solution'],false); 0202 warned(1)=true; 0203 end 0204 else 0205 if warned(2)==false 0206 dispEM(['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 random reactions in the solution'],false); 0207 warned(2)=true; 0208 end 0209 end 0210 I=find(abs(solution)>10^-8); 0211 rxnToRemove=I(randsample(numel(I),1)); 0212 end 0213 end 0214 end 0215 end 0216 removedRxns=[removedRxns;model.rxns(rxnToRemove)]; 0217 if printReport 0218 fprintf(['\tRemoved: ' model.rxns{rxnToRemove} '\n']); 0219 end 0220 model=removeRxns(model,rxnToRemove); 0221 balanceStructure.balanceStatus(rxnToRemove)=[]; 0222 balanceStructure.leftComp(rxnToRemove,:)=[]; 0223 balanceStructure.rightComp(rxnToRemove,:)=[]; 0224 end 0225 end 0226 0227 newModel=model; 0228 end