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 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

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 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

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005