removeBadRxns Iteratively removes reactions which enable production of some metabolite without any uptake model a model structure. For the intented function, the model shouldn't allow for any uptake. 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 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) 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 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 of metabolites without any consumption. 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 anything from nothing. This is done by repeatedly calling makeSomething, checking if any of the involved reactions are elementally unbalanced, remove one of them, and then iterating until no metabolites can be synthesized. Usage: [newModel removedRxns]=removeBadRxns(model,rxnRules,... ignoreMets,isNames,refModel) Rasmus Agren, 2012-03-20
0001 function [newModel removedRxns]=removeBadRxns(model,rxnRules,ignoreMets,isNames,balanceElements,refModel) 0002 % removeBadRxns 0003 % Iteratively removes reactions which enable production of some 0004 % metabolite without any uptake 0005 % 0006 % model a model structure. For the intented function, 0007 % the model shouldn't allow for any uptake. The 0008 % 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 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 % 0039 % newModel a model structure after the problematic 0040 % reactions have been deleted 0041 % removedRxns a cell array with the reactions that were 0042 % removed 0043 % 0044 % The purpose of this function is to remove reactions which enable 0045 % production of metabolites even when exchange reactions aren't used. 0046 % Many models, especially if they are automatically inferred from 0047 % databases, will have unbalanced reactions which allow for 0048 % net-production of metabolites without any consumption. A common reason 0049 % for this is when general compounds have different meaning in different 0050 % reactions (as DNA has in these two reactions). 0051 % dATP + dGTP + dCTP + dTTP <=> DNA + 4 PPi 0052 % 0.25 dATP + 0.25 dGTP + 0.25 dCTP + 0.25 dTTP <=> DNA + PPi 0053 % Reactions that are problematic like this are always elementally 0054 % unbalanced, but it is not always that you would like to exclude all 0055 % unbalanced reactions from your model. 0056 % This function tries to remove as few problematic reactions as possible 0057 % so that the model cannot produce anything from nothing. This is done by 0058 % repeatedly calling makeSomething, checking if any of the involved 0059 % reactions are elementally unbalanced, remove one of them, and then 0060 % iterating until no metabolites can be synthesized. 0061 % 0062 % Usage: [newModel removedRxns]=removeBadRxns(model,rxnRules,... 0063 % ignoreMets,isNames,refModel) 0064 % 0065 % Rasmus Agren, 2012-03-20 0066 % 0067 0068 if nargin<2 0069 rxnRules=1; 0070 end 0071 if nargin<3 0072 ignoreMets=[]; 0073 end 0074 if nargin<4 0075 isNames=false; 0076 end 0077 if nargin<5 0078 balanceElements={'C';'P';'S';'N';'O'}; 0079 end 0080 if nargin<6 0081 refModel=[]; 0082 end 0083 0084 %Check if the model has open exchange reactions and print a warning in that case 0085 if ~isfield(model,'unconstrained') 0086 [crap I]=getExchangeRxns(model); 0087 if any(I) 0088 if any(model.lb(I)~=0) || any(model.ub(I)~=0) 0089 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'); 0090 end 0091 end 0092 end 0093 0094 %Check that the model is feasible 0095 sol=solveLP(model); 0096 if isempty(sol.f) 0097 throw(MException('','The model is not feasible. Consider removing lower bounds (such as ATP maintenance)')); 0098 end 0099 0100 %Check that the reference model is feasible 0101 if any(refModel) 0102 sol=solveLP(refModel); 0103 if isempty(sol.f) 0104 throw(MException('','The reference model is not feasible.')); 0105 end 0106 end 0107 0108 %Initialize some stuff 0109 removedRxns={}; 0110 balanceStructure=getElementalBalance(model); 0111 0112 %Check which elements to balance for 0113 if ~isempty(setdiff(upper(balanceElements),{'C','P','S','N','R','O','H'})) 0114 throw(MException('','Could not recognize all elements to balance for')); 0115 end 0116 bal=ismember({'C','P','S','N','R','O','H'},upper(balanceElements)); 0117 0118 %Main loop 0119 while 1 0120 %Make some metabolite using as few reactions as possible 0121 [solution metabolite]=makeSomething(model,ignoreMets,isNames); 0122 0123 %If no soulution could be found, then finish 0124 if isempty(solution) 0125 break; 0126 end 0127 0128 %Find all reactions that are unbalanced and still carry flux 0129 I=find(abs(solution)>10^-8 & balanceStructure.balanceStatus>=0 & ~all(balanceStructure.leftComp(:,bal)==balanceStructure.rightComp(:,bal),2)); 0130 0131 %If there are unbalanced rxns then delete one of them and iterate 0132 if any(I) 0133 rxnToRemove=I(randsample(numel(I),1)); 0134 else 0135 %If there are no unbalanced rxns in the solution 0136 if rxnRules==1 0137 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']); 0138 break; 0139 else 0140 %Find reactions which are not checked for mass balancing, but 0141 %that still carry flux 0142 I=find(abs(solution)>10^-8 & balanceStructure.balanceStatus<0); 0143 0144 %If there are any such reactions, remove one of them and 0145 %iterate 0146 if any(I) 0147 rxnToRemove=I(randsample(numel(I),1)); 0148 else 0149 if rxnRules==2 0150 %This happens when all reactions used are balanced 0151 %according to the metabolite formulas. This cannot be, and 0152 %indicates that one or more of the formulas are wrong. 0153 %Print a warning and delete any reaction with flux 0154 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']); 0155 break; 0156 else 0157 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']); 0158 I=find(abs(solution)>10^-8); 0159 rxnToRemove=I(randsample(numel(I),1)); 0160 end 0161 end 0162 end 0163 end 0164 removedRxns=[removedRxns;model.rxns(rxnToRemove)]; 0165 model=removeRxns(model,rxnToRemove); 0166 balanceStructure.balanceStatus(rxnToRemove)=[]; 0167 balanceStructure.leftComp(rxnToRemove,:)=[]; 0168 balanceStructure.rightComp(rxnToRemove,:)=[]; 0169 end 0170 0171 newModel=model; 0172 end