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)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 if nargin<2 0073 rxnRules=1; 0074 end 0075 if nargin<3 0076 ignoreMets=[]; 0077 elseif ~islogical(ignoreMets) && ~isnumeric(ignoreMets) 0078 ignoreMets=convertCharArray(ignoreMets); 0079 end 0080 if nargin<4 0081 isNames=false; 0082 end 0083 if nargin<5 0084 balanceElements={'C';'P';'S';'N';'O'}; 0085 else 0086 balanceElements=convertCharArray(balanceElements); 0087 end 0088 if nargin<6 0089 refModel=[]; 0090 else 0091 if ~isempty(refModel) 0092 EM='The feature to supply a reference model is currently not supported'; 0093 dispEM(EM,false); 0094 end 0095 end 0096 if nargin<7 0097 ignoreIntBounds=false; 0098 end 0099 if nargin<8 0100 printReport=false; 0101 end 0102 0103 %Check if the model has open exchange reactions and print a warning in that 0104 %case 0105 if ~isfield(model,'unconstrained') 0106 [~, I]=getExchangeRxns(model); 0107 if any(I) 0108 if any(model.lb(I)~=0) || any(model.ub(I)~=0) 0109 EM='The model contains open exchange reactions. This is not the intended use of this function. Consider importing your model using importModel(filename,false)'; 0110 dispEM(EM,false); 0111 end 0112 end 0113 end 0114 0115 %Check that the model is feasible 0116 sol=solveLP(model); 0117 if isempty(sol.f) 0118 EM='The model is not feasible. Consider removing lower bounds (such as ATP maintenance)'; 0119 dispEM(EM); 0120 end 0121 0122 %Check that the reference model is feasible 0123 if any(refModel) 0124 sol=solveLP(refModel); 0125 if isempty(sol.f) 0126 EM='The reference model is not feasible'; 0127 dispEM(EM); 0128 end 0129 end 0130 0131 %Initialize some stuff 0132 removedRxns={}; 0133 balanceStructure=getElementalBalance(model); 0134 0135 %Check which elements to balance for 0136 if ~isempty(setdiff(balanceElements,balanceStructure.elements.abbrevs)) 0137 EM='Could not recognize all elements to balance for'; 0138 dispEM(EM); 0139 end 0140 bal=ismember(balanceStructure.elements.abbrevs,balanceElements); 0141 0142 %Main loop. First run for makeSomething, second for consumeSomething 0143 warned=false(2,1); %This is to prevent the same warning being printed multiple times if rxnRules==3 0144 for i=1:2 0145 while 1 0146 %Make some metabolite using as few reactions as possible 0147 if i==1 0148 [solution, metabolite]=makeSomething(model,ignoreMets,isNames,false,true,[],ignoreIntBounds); 0149 if ~isempty(solution) 0150 if printReport 0151 fprintf(['Can make: ' model.metNames{metabolite(1)} '\n']); 0152 end 0153 else 0154 %If no solution could be found, then finish 0155 break; 0156 end 0157 else 0158 [solution, metabolite]=consumeSomething(model,ignoreMets,isNames,false,[],ignoreIntBounds); 0159 if ~isempty(solution) 0160 if printReport 0161 fprintf(['Can consume: ' model.metNames{metabolite(1)} '\n']); 0162 end 0163 else 0164 %If no solution could be found, then finish 0165 break; 0166 end 0167 end 0168 0169 %Find all reactions that are unbalanced and still carry flux 0170 I=find(abs(solution)>10^-8 & balanceStructure.balanceStatus>=0 & ~all(balanceStructure.leftComp(:,bal)==balanceStructure.rightComp(:,bal),2)); 0171 0172 %If there are unbalanced rxns then delete one of them and iterate 0173 if any(I) 0174 rxnToRemove=I(randperm(numel(I),1)); 0175 else 0176 %If there are no unbalanced rxns in the solution 0177 if rxnRules==1 0178 if i==1 0179 EM=['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']; 0180 dispEM(EM,false); 0181 else 0182 EM=['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']; 0183 dispEM(EM,false); 0184 end 0185 break; 0186 else 0187 %Find reactions which are not checked for mass balancing, 0188 %but that still carry flux 0189 I=find(abs(solution)>10^-8 & balanceStructure.balanceStatus<0); 0190 0191 %If there are any such reactions, remove one of them and 0192 %iterate 0193 if any(I) 0194 rxnToRemove=I(randperm(numel(I),1)); 0195 else 0196 if rxnRules==2 0197 %This happens when all reactions used are balanced 0198 %according to the metabolite formulas. This cannot 0199 %be, and indicates that one or more of the formulas 0200 %are wrong. Print a warning and delete any reaction 0201 %with flux 0202 if i==1 0203 EM=['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']; 0204 dispEM(EM,false); 0205 else 0206 EM=['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']; 0207 dispEM(EM,false); 0208 end 0209 break; 0210 else 0211 if i==1 0212 if warned(1)==false 0213 EM=['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']; 0214 dispEM(EM,false); 0215 warned(1)=true; 0216 end 0217 else 0218 if warned(2)==false 0219 EM=['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']; 0220 dispEM(EM,false); 0221 warned(2)=true; 0222 end 0223 end 0224 I=find(abs(solution)>10^-8); 0225 rxnToRemove=I(randperm(numel(I),1)); 0226 end 0227 end 0228 end 0229 end 0230 removedRxns=[removedRxns;model.rxns(rxnToRemove)]; 0231 if printReport 0232 fprintf(['\tRemoved: ' model.rxns{rxnToRemove} '\n']); 0233 end 0234 model=removeReactions(model,rxnToRemove); 0235 balanceStructure.balanceStatus(rxnToRemove)=[]; 0236 balanceStructure.leftComp(rxnToRemove,:)=[]; 0237 balanceStructure.rightComp(rxnToRemove,:)=[]; 0238 end 0239 end 0240 0241 newModel=model; 0242 end