Home > RAVEN > simplifyModel.m

simplifyModel

PURPOSE ^

simplifyModel

SYNOPSIS ^

function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, reservedRxns, suppressWarnings)

DESCRIPTION ^

 simplifyModel
   Simplifies a model by deleting reactions/metabolites

   model                 a model structure
   deleteUnconstrained   delete metabolites marked as unconstrained (opt, default true)
   deleteDuplicates      delete all but one of duplicate reactions (opt, default false)
   deleteZeroInterval    delete reactions that are constrained to zero flux (opt, default false)
   deleteInaccessible    delete dead end reactions (opt, default false)
   deleteMinMax          delete reactions that cannot carry a flux by trying
                         to minimize/maximize the flux through that
                         reaction. May be time consuming (opt, default false)
   groupLinear           group linear pathways (opt, default false)
   reservedRxns          cell array with reaction IDs that are not allowed to be
                         removed (opt)
   suppressWarnings      true if warnings should be suppressed (opt,
                         default false)

   reducedModel          an updated model structure
   deletedReactions      a cell array with the IDs of all deleted reactions
   deletedMetabolites    a cell array with the IDs of all deleted
                         metabolites

   This function is for reducing the model size by removing
   reactions and associated metabolites that cannot carry flux. It can also
   be used for identifying different types of gaps.

   Usage: [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,...
           deleteUnconstrained, deleteDuplicates, deleteZeroInterval,...
           deleteInaccessible, deleteMinMax, groupLinear, reservedRxns,...
           suppressWarnings)

   Rasmus Agren, 2013-06-21

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,...
0002     deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, reservedRxns, suppressWarnings)
0003 % simplifyModel
0004 %   Simplifies a model by deleting reactions/metabolites
0005 %
0006 %   model                 a model structure
0007 %   deleteUnconstrained   delete metabolites marked as unconstrained (opt, default true)
0008 %   deleteDuplicates      delete all but one of duplicate reactions (opt, default false)
0009 %   deleteZeroInterval    delete reactions that are constrained to zero flux (opt, default false)
0010 %   deleteInaccessible    delete dead end reactions (opt, default false)
0011 %   deleteMinMax          delete reactions that cannot carry a flux by trying
0012 %                         to minimize/maximize the flux through that
0013 %                         reaction. May be time consuming (opt, default false)
0014 %   groupLinear           group linear pathways (opt, default false)
0015 %   reservedRxns          cell array with reaction IDs that are not allowed to be
0016 %                         removed (opt)
0017 %   suppressWarnings      true if warnings should be suppressed (opt,
0018 %                         default false)
0019 %
0020 %   reducedModel          an updated model structure
0021 %   deletedReactions      a cell array with the IDs of all deleted reactions
0022 %   deletedMetabolites    a cell array with the IDs of all deleted
0023 %                         metabolites
0024 %
0025 %   This function is for reducing the model size by removing
0026 %   reactions and associated metabolites that cannot carry flux. It can also
0027 %   be used for identifying different types of gaps.
0028 %
0029 %   Usage: [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,...
0030 %           deleteUnconstrained, deleteDuplicates, deleteZeroInterval,...
0031 %           deleteInaccessible, deleteMinMax, groupLinear, reservedRxns,...
0032 %           suppressWarnings)
0033 %
0034 %   Rasmus Agren, 2013-06-21
0035 %
0036 
0037 if nargin<2
0038     deleteUnconstrained=true;
0039 end
0040 if nargin<3
0041     deleteDuplicates=false;
0042 end
0043 if nargin<4
0044     deleteZeroInterval=false;
0045 end
0046 if nargin<5
0047     deleteInaccessible=false;
0048 end
0049 if nargin<6
0050     deleteMinMax=false;
0051 end
0052 if nargin<7
0053     groupLinear=false;
0054 end
0055 if nargin<8
0056     reservedRxns=[];
0057 end
0058 if nargin<9
0059     suppressWarnings=false;
0060 end
0061 
0062 reducedModel=model;
0063 deletedReactions={};
0064 deletedMetabolites={};
0065 
0066 if deleteUnconstrained==true
0067     if isfield(reducedModel,'unconstrained')
0068         %Remove unbalanced metabolites
0069         deletedMetabolites=reducedModel.mets(reducedModel.unconstrained~=0);
0070         reducedModel=removeMets(reducedModel,reducedModel.unconstrained~=0);
0071         reducedModel=rmfield(reducedModel,'unconstrained');
0072     end
0073 end
0074 
0075 if deleteDuplicates==true
0076     %Delete all but the last occurrence of duplicate reactions. The
0077     %reactions must have the same bounds, reversibility, and objective
0078     %coefficient to be regarded as duplicate
0079     [reducedModel rxnsToDelete]=contractModel(reducedModel);
0080     deletedReactions=[deletedReactions; rxnsToDelete];
0081 end
0082 
0083 if deleteZeroInterval==true
0084     %Find all reactions where both LB and UB are 0
0085     zeroIntervalReactions=and(reducedModel.lb==0,reducedModel.ub==0);
0086     
0087     rxnsToDelete=setdiff(reducedModel.rxns(zeroIntervalReactions),reservedRxns);
0088     deletedReactions=[deletedReactions; rxnsToDelete];   
0089     
0090     %Remove reactions
0091     reducedModel=removeRxns(reducedModel,rxnsToDelete);
0092     
0093     %Find metabolites that no longer are used and delete them
0094     notInUse=sum(reducedModel.S~=0,2)==0;
0095     deletedMetabolites=[deletedMetabolites;reducedModel.mets(notInUse)];
0096     
0097     %Remove metabolites
0098     reducedModel=removeMets(reducedModel,notInUse);
0099 end
0100 
0101 if deleteInaccessible==true
0102     %Print a warning if exchange metabolites haven't been deleted yet. This
0103     %often means that the only allowed products are the metabolites that
0104     %are taken up be the system
0105     if isfield(reducedModel,'unconstrained') && suppressWarnings==false
0106        fprintf('WARNING: Removing dead-end reactions before removing exchange metabolites\n'); 
0107     end
0108     
0109     while true
0110         %Find all metabolites that are only reactants or only products.
0111         %Delete those metabolites and reactions. This is done until no more
0112         %such reactions are found
0113         
0114         %Construct a stoichiometric matrix where all reactions are written
0115         %as reversible
0116         
0117         %Add fake exchange reactions if model.b~=0
0118         in=any(reducedModel.b<0,2);
0119         out=any(reducedModel.b>0,2);
0120         I=speye(numel(reducedModel.mets));
0121         revS=[reducedModel.S,reducedModel.S(:,reducedModel.rev~=0)*-1 I(:,in) I(:,out)*-1];
0122         
0123         metUsage=sum(abs(revS')>0);
0124         onlyProducts=sum(revS'>0) == metUsage;
0125         onlyReactants=sum(revS'<0) == metUsage;
0126         
0127         %Also remove metabolites that only participate in one reversible reaction
0128         %Don't remove the ones in one reversible reaction if is also has a
0129         %non-zero coefficient in model.b
0130         notInUse=onlyProducts | onlyReactants | (sum(abs(reducedModel.S')>0)<=1 & (~in & ~out)');
0131         deletedRxn=false;
0132         if any(notInUse)        
0133             %Find their corresponding reactions
0134             rxnsNotInUse=sum(abs(reducedModel.S(notInUse,:))>0,1)>0;
0135             rxnsToDelete=setdiff(reducedModel.rxns(rxnsNotInUse),reservedRxns);
0136             deletedReactions=[deletedReactions; rxnsToDelete];
0137             
0138             %Remove reactions
0139             reducedModel=removeRxns(reducedModel,rxnsToDelete);
0140 
0141             %Remove metabolites. Recalculate since it could be that some
0142             %cannot be deleted due to reserved rxns
0143             notInUse=sum(reducedModel.S~=0,2)==0;
0144             deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0145             reducedModel=removeMets(reducedModel,notInUse);
0146             
0147             %It could happen that it just finds reserved reactions and gets
0148             %stuck
0149             if ~isempty(rxnsToDelete)
0150                 deletedRxn=true;
0151             end
0152         else
0153             break;
0154         end
0155         if deletedRxn==false
0156             break;
0157         end
0158     end
0159 end
0160 
0161 if deleteMinMax==true
0162     %Get reactions that can't carry fluxes. This should be done
0163     %algebraically if possible
0164     I=~haveFlux(reducedModel);
0165     
0166     %Remove reactions
0167     rxnsToDelete=setdiff(reducedModel.rxns(I),reservedRxns);
0168     deletedReactions=[deletedReactions; rxnsToDelete];
0169     reducedModel=removeRxns(reducedModel,rxnsToDelete);
0170             
0171     %Remove metabolites
0172     notInUse=sum(reducedModel.S~=0,2)==0;
0173     deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0174     reducedModel=removeMets(reducedModel,notInUse);
0175 end
0176 
0177 %Checks that all reactions are irreversible. Might be fixed in the future.
0178 if groupLinear==true
0179     if ~any(reducedModel.rev)
0180         if  suppressWarnings==false
0181             fprintf('NOTE: You have chosen to group linear reactions. This option does not keep track of gene/reaction associations when reactions are merged. Deleting all gene information...\n');
0182         end
0183         bannedIndexes=getIndexes(reducedModel,reservedRxns,'rxns');
0184         while 1
0185             %Select all metabolites that are only present as reactants/products
0186             %in one reaction
0187             singleNegative=find(sum(reducedModel.S'<0)==1);
0188             singlePositive=find(sum(reducedModel.S'>0)==1);
0189    
0190             %Retrieve the common metabolites
0191             common=intersect(singleNegative,singlePositive);
0192             
0193             %Merge the reactions for each of these metabolites
0194             mergedSomeRxn=false;
0195             for i=1:numel(common)
0196                 involvedRxns=find(reducedModel.S(common(i),:));
0197                 
0198                 %Check so that it doesn't try to merge reactions that are
0199                 %restricted. The second check is needed because it could be
0200                 %that some of the metabolites have already been merged in
0201                 %this iteration so that the reactions no longer exist
0202                 if isempty(intersect(bannedIndexes,involvedRxns)) && any(involvedRxns)
0203                     %This is the ratio that described how many times the second
0204                     %reaction should be multiplied before being merged with
0205                     %the first
0206                     stoichRatio=abs(reducedModel.S(common(i),involvedRxns(1))/reducedModel.S(common(i),involvedRxns(2)));
0207                     reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio;
0208                     reducedModel.S(common(i),involvedRxns(1))=0;
0209                     
0210                     %Change the name and id to reflect the merging
0211                     reducedModel.rxns{involvedRxns(1)}=[reducedModel.rxns{involvedRxns(1)} '_' reducedModel.rxns{involvedRxns(2)}];
0212                     reducedModel.rxnNames{involvedRxns(1)}=reducedModel.rxns{involvedRxns(1)};
0213                     
0214                     %Remove the second reaction by setting all coefficients
0215                     %to 0
0216                     reducedModel.S(:,involvedRxns(2))=0;
0217                     mergedSomeRxn=true;
0218                 end
0219             end
0220             
0221             %If it couldn't merge anything then exit the loop
0222             if mergedSomeRxn==false
0223                 break;
0224             end
0225         end
0226         
0227         %Now delete all reactions that involve no metabolites
0228         I=sum(reducedModel.S~=0)==0;
0229         
0230         %Remove reactions
0231         deletedReactions=[deletedReactions; reducedModel.rxns(I)];
0232         reducedModel=removeRxns(reducedModel,find(I));
0233             
0234         %Remove metabolites
0235         notInUse=sum(reducedModel.S~=0,2)==0;
0236         deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0237         reducedModel=removeMets(reducedModel,notInUse);
0238         
0239         reducedModel.genes={};
0240         reducedModel.rxnGeneMat=sparse(numel(reducedModel.rxns),0);
0241         reducedModel.grRules(:)={''};
0242         
0243         if isfield(reducedModel,'geneShortNames')
0244             reducedModel.geneShortNames={};
0245         end
0246         if isfield(reducedModel,'geneMiriams')
0247             reducedModel.geneMiriams={};
0248         end
0249     else
0250         throw(MException('','Cannot group reactions for a model that has reversible reactions'));
0251     end
0252 end
0253 
0254 end

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