Home > RAVEN > removeBadRxns.m

removeBadRxns

PURPOSE ^

removeBadRxns

SYNOPSIS ^

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

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005