Home > RAVEN > removeBadRxns.m

removeBadRxns

PURPOSE ^

removeBadRxns

SYNOPSIS ^

function [newModel removedRxns]=removeBadRxns(model,rxnRules,ignoreMets,isNames,balanceElements,refModel,ignoreIntBounds,printReport)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005