contractModel

 contractModel
   Contracts a model by grouping all identical reactions. Similar to the
   deleteDuplicates part in simplifyModel but more care is taken here
   when it comes to gene associations. If the duplicated reactions have
   '_EXP_*' suffixes (where * is a digit), then the model is assumed to
   have been passed through expandModel, and these suffixes are removed
   here.

   model                  a model structure
   distReverse            distinguish reactions with same metabolites
                          but different reversibility as different
                          reactions (opt, default true)
   mets                   string or cell array of strings with metabolite
                          identifiers, whose involved reactions should be
                          checked for duplication (opt, by default all
                          reactions are considered) (option is used by
                          replaceMets)

   reducedModel           a model structure without duplicate reactions
   removedRxns            cell array for the removed duplicate reactions
   indexedDuplicateRxns   indexed cell array for the removed duplicate
                          reactions (multiple valuess separated by semicolon)

   NOTE: This code might not work for advanced grRules strings
         that involve nested expressions of 'and' and 'or'.

   Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse)

0001 function [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse,mets)
0002 % contractModel
0003 %   Contracts a model by grouping all identical reactions. Similar to the
0004 %   deleteDuplicates part in simplifyModel but more care is taken here
0005 %   when it comes to gene associations. If the duplicated reactions have
0006 %   '_EXP_*' suffixes (where * is a digit), then the model is assumed to
0007 %   have been passed through expandModel, and these suffixes are removed
0008 %   here.
0009 %
0010 %   model                  a model structure
0011 %   distReverse            distinguish reactions with same metabolites
0012 %                          but different reversibility as different
0013 %                          reactions (opt, default true)
0014 %   mets                   string or cell array of strings with metabolite
0015 %                          identifiers, whose involved reactions should be
0016 %                          checked for duplication (opt, by default all
0017 %                          reactions are considered) (option is used by
0018 %                          replaceMets)
0019 %
0020 %   reducedModel           a model structure without duplicate reactions
0021 %   removedRxns            cell array for the removed duplicate reactions
0022 %   indexedDuplicateRxns   indexed cell array for the removed duplicate
0023 %                          reactions (multiple valuess separated by semicolon)
0024 %
0025 %   NOTE: This code might not work for advanced grRules strings
0026 %         that involve nested expressions of 'and' and 'or'.
0027 %
0028 %   Usage: [reducedModel, removedRxns, indexedDuplicateRxns]=contractModel(model,distReverse)
0029 
0030 if nargin<2
0031     distReverse=true;
0032 end
0033 
0034 %First sort the model so that reversible reactions are in the same
0035 %direction
0036 modelS=sortModel(model);
0037 
0038 %Get a list of duplicate reactions
0039 if distReverse
0040     x=[modelS.S; model.rev']';
0041 else
0042     x=modelS.S';
0043 end
0044 if nargin>2
0045     toCheck = getIndexes(model,mets,'mets');
0046     [~, toCheck] = find(modelS.S(toCheck,:));
0047     x = [x, (1:size(x,1))']; % Make all other rxns unique by additional column
0048     x(toCheck,end) = 0;
0049 end
0050 
0051 [~,I,J] = unique(x,'rows','first');
0052 
0053 %Initialize cell array of indexedDuplicateRxns
0054 indexedDuplicateRxns=cell(numel(model.rxns),1);
0055 indexedDuplicateRxns(:)={''};
0056 
0057 duplicateRxns=setdiff(1:numel(model.rxns),I);
0058 mergeTo=I(J(duplicateRxns));
0059 
0060 mergedRxns=unique(mergeTo);
0061 
0062 %Now add all the info from this one. Print a warning if they have different
0063 %bounds or objective function coefficients. Uses the widest bounds and
0064 %largest magnitude of objective coefficient
0065 for i=1:numel(mergedRxns)
0066     duplRxn=transpose([mergedRxns(i),duplicateRxns(mergeTo==mergedRxns(i))]);
0067     if numel(unique(model.lb(duplRxn)))>1
0068         EM=['Duplicates of reaction ' model.rxns{mergedRxns(i)} ' have different lower bound. Uses the most negative/smallest lower bound'];
0069         dispEM(EM,false);
0070         model.lb(mergedRxns(i))=min(model.lb(duplRxn));
0071     end
0072     if numel(unique(model.ub(duplRxn)))>1
0073         EM=['Duplicates of reaction ' model.rxns{mergedRxns(i)} ' have different upper bound. Uses the most positive/largest upper bound'];
0074         dispEM(EM,false);
0075         model.ub(mergedRxns(i))=max(model.ub(duplRxn));
0076     end
0077     if numel(unique(model.c(duplRxn)))>1
0078         EM=['Duplicates of reaction ' model.rxns{mergedRxns(i)} ' has a different objective function coefficient. Uses the largest coefficient'];
0079         dispEM(EM,false);
0080         model.c(mergedRxns(i))=max(model.c(duplRxn));
0081     end
0082     if isfield(model,'grRules') && any(~isempty(model.grRules(duplRxn)))
0083         rules=model.grRules(duplRxn);
0084         allRules={};
0085         for j=1:numel(rules)
0086             rules{j}=ignoreORinComplex(rules{j});
0087             allRules=[allRules regexp(rules{j},' or ','split')];
0088         end
0089         allRules=unique(allRules);
0090         allRules=strrep(allRules,'__OR__',' or ');
0091         andRules=contains(allRules,' and ');
0092         allRules(andRules)=strcat('(',allRules(andRules),')');
0093         if numel(allRules)==1
0094             model.grRules{mergedRxns(i)}=allRules{1};
0095         else
0096             model.grRules{mergedRxns(i)}=strjoin(allRules,' or ');
0097         end
0098     end    
0099     if isfield(model,'eccodes') && any(~isempty(model.eccodes(duplRxn)))
0100         codes=model.eccodes(duplRxn);
0101         allCodes={};
0102         for j=1:numel(codes)
0103             allCodes=[allCodes regexp(codes{j},';','split')];
0104         end
0105         allCodes=unique(allCodes);
0106         if numel(allCodes)==1
0107             model.eccodes{mergedRxns(i)}=allCodes{1};
0108         else
0109             model.eccodes{mergedRxns(i)}=strjoin(allCodes,';');
0110         end
0111     end
0112     %Generate indexedDuplicateRxns cell array
0113     if numel(duplRxn)==2
0114         indexedDuplicateRxns{duplRxn(1)}=model.rxns{duplRxn(2)};
0115     else
0116         indexedDuplicateRxns{duplRxn(1)}=strjoin(model.rxns(duplRxn(2:end)),';');
0117     end
0118     %If all reactions have _EXP_* (where * is digit) as suffix, then the
0119     %model had passed through expandModel(), and these suffixes are removed.
0120     if all(~cellfun(@isempty,regexp(model.rxns(duplRxn),'_EXP_\d+$')))
0121         model.rxns(duplRxn(1))=regexprep(model.rxns(duplRxn(1)),'_EXP_\d+$','');
0122     end
0123 end
0124 
0125 %Delete the duplicate reactions
0126 reducedModel=removeReactions(model,duplicateRxns);
0127 removedRxns=model.rxns(duplicateRxns);
0128 [~, index]=ismember(reducedModel.rxns,model.rxns);
0129 indexedDuplicateRxns=indexedDuplicateRxns(index);
0130 
0131 if isfield(reducedModel,'rxnGeneMat')
0132     %Fix grRules and reconstruct rxnGeneMat
0133     [grRules,rxnGeneMat] = standardizeGrRules(reducedModel,true);
0134     reducedModel.grRules = grRules;
0135     reducedModel.rxnGeneMat = rxnGeneMat;
0136 end
0137 end
0138 
0139 function grRule = ignoreORinComplex(grRule)
0140 %In a grRule, if OR relationship is nested in an AND relationship, then
0141 %obfuscate the OR before splitting isoenzymes
0142 grRule=['(' grRule ')'];
0143 brOpen=strfind(grRule,'(');
0144 brClose=strfind(grRule,')');
0145 andPos=strfind(grRule,' and ');
0146 %Find opening bracket closest before AND
0147 stillCapturing = 0;
0148 for i=1:numel(andPos)
0149     searchPos = andPos(i);
0150     while stillCapturing == 0
0151         closestOpen = brOpen(max(find(brOpen<searchPos)));
0152         inbetweenClose = brClose(brClose<searchPos & brClose>closestOpen);
0153         if ~isempty(inbetweenClose)
0154             searchPos=max(inbetweenClose);
0155         else
0156             stillCapturing = 1;
0157             beginPos = closestOpen;
0158         end
0159     end
0160     stillCapturing = 0;
0161     searchPos = andPos(i);
0162     while stillCapturing == 0
0163         closestClose = brClose(min(find(brClose>searchPos)));
0164         inbetweenOpen = brOpen(brOpen>searchPos & brOpen<closestOpen);
0165         if ~isempty(inbetweenOpen)
0166             searchPos=min(closestClose);
0167         else
0168             stillCapturing = 1;
0169             endPos = closestClose;
0170         end
0171     end
0172     replacePart=regexprep(grRule(beginPos:endPos),' or ','__OR__');
0173     grRule=strrep(grRule,grRule(beginPos:endPos),replacePart);
0174 end
0175 grRule=grRule(2:end-1);
0176 end