mergeLinear

 mergeLinear
   Simplifies a model by merging rxns with linear dependencies. There are two
   differences from the function in simplifyModel:
   1) Here we return which reactions were merged
   2) Here we allow reversible reactions to be merged as well, which will 
      probably reduce the size of the model more.

   model                 a model structure
   noMergeRxns           Cell array with reaction IDs that are not allowed to be
                         merged

   reducedModel          an updated model structure
   origRxnIds            Cell array: original rxn ids, used together with the groupIds variable
   groupIds              Vector of ids: merged rxns have the same id, same length as origRxnIds
                         0 means a reaction was not merged.
   reversedRxns          Vector of booleans, saying for each rxn if it is reversed in the combined rxns.

 Usage: [reducedModel,origRxnIds,groupIds]=mergeLinear(model,noMergeRxns)

0001 function [reducedModel,origRxnIds,groupIds,reversedRxns]=mergeLinear(model, noMergeRxns)
0002 % mergeLinear
0003 %   Simplifies a model by merging rxns with linear dependencies. There are two
0004 %   differences from the function in simplifyModel:
0005 %   1) Here we return which reactions were merged
0006 %   2) Here we allow reversible reactions to be merged as well, which will
0007 %      probably reduce the size of the model more.
0008 %
0009 %   model                 a model structure
0010 %   noMergeRxns           Cell array with reaction IDs that are not allowed to be
0011 %                         merged
0012 %
0013 %   reducedModel          an updated model structure
0014 %   origRxnIds            Cell array: original rxn ids, used together with the groupIds variable
0015 %   groupIds              Vector of ids: merged rxns have the same id, same length as origRxnIds
0016 %                         0 means a reaction was not merged.
0017 %   reversedRxns          Vector of booleans, saying for each rxn if it is reversed in the combined rxns.
0018 %
0019 % Usage: [reducedModel,origRxnIds,groupIds]=mergeLinear(model,noMergeRxns)
0020 
0021 
0022 reducedModel=model;
0023 
0024 
0025 reducedModel.genes={};
0026 reducedModel.rxnGeneMat=sparse(numel(reducedModel.rxns),0);
0027 reducedModel.grRules(:)={''};
0028 
0029 if isfield(reducedModel,'geneShortNames')
0030     reducedModel.geneShortNames={};
0031 end
0032 if isfield(reducedModel,'proteins')
0033     reducedModel.proteins={};
0034 end
0035 if isfield(reducedModel,'geneMiriams')
0036     reducedModel.geneMiriams={};
0037 end
0038 if isfield(reducedModel,'geneComps')
0039     reducedModel.geneComps=[];
0040 end
0041 
0042 nextGroupId = 1;
0043 origRxnIds = reducedModel.rxns;
0044 groupIds = zeros(numel(reducedModel.rxns),1);
0045 reversedRxns = false(numel(reducedModel.rxns),1);
0046 
0047 %Loop through and iteratively group linear reactions
0048 while 1
0049     %Get the banned reaction indexes. Note that the indexes will change
0050     %in each iteration, but the names will not as they won't be merged
0051     %with any other reaction
0052     bannedIndexes=getIndexes(reducedModel,noMergeRxns,'rxns');
0053 
0054     %Select all metabolites that are only present as reactants/products
0055     %in one reaction
0056     twoNonZero = find(sum(reducedModel.S ~= 0, 2) == 2);
0057 
0058     mergedSome=false;
0059 
0060     %Loop through each of them and see if the reactions should be
0061     %merged
0062     for i=1:numel(twoNonZero)
0063         involvedRxns=find(reducedModel.S(twoNonZero(i),:));
0064         %Check that we can have one positive and one negative
0065         %rev is compared rather than used as-is: the other operand is a
0066         %sparse logical, and MATLAB refuses to combine one with an integer
0067         %array ("Sparse integer array arithmetic operations are not
0068         %supported"). A model built outside RAVEN can carry an integer rev.
0069         pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0);
0070         neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0);
0071         
0072                 
0073         %Check so that one or both of the reactions haven't been merged
0074         %already
0075         if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1
0076             %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible,
0077             %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise
0078             if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ...
0079                (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0)
0080                 involvedRxns = flip(involvedRxns);
0081             end
0082             
0083             %first make sure the first reaction is producing the metabolite
0084             if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0
0085                 %it is not producing the metabolite - fix that
0086                 %first choice: use the second reaction as producer if it is producing
0087                 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0
0088                     involvedRxns = flip(involvedRxns);%make the second the first
0089                 else
0090                     %now we know that the second reaction is not producing, so we can safely try to make the first a producer
0091                     if reducedModel.rev(involvedRxns(1)) == 1
0092                         [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns);
0093                     else %ok, finally try to flip the second reaction
0094                         if reducedModel.rev(involvedRxns(2)) == 1
0095                             [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns);
0096                             involvedRxns = flip(involvedRxns);%make the second the first
0097                         else
0098                             error('We should never end up here!');
0099                         end
0100                     end
0101                 end
0102             end
0103             %Now, make sure the second rxn is a consumer
0104             if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0
0105                 if reducedModel.rev(involvedRxns(2)) == 1
0106                     [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns);
0107                 else
0108                     error('We should never end up here!');
0109                 end
0110             end
0111 
0112             %Calculate how many times the second reaction has to be
0113             %multiplied before being merged with the first
0114             stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2)));
0115 
0116             %Add the second to the first
0117             reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio;
0118 
0119             %Clear the second reaction
0120             reducedModel.S(:,involvedRxns(2))=0;
0121 
0122             %This is to prevent numerical issues. It should be 0
0123             %already
0124             reducedModel.S(twoNonZero(i),involvedRxns(1))=0;
0125 
0126             %At this point the second reaction is certain to be deleted
0127             %in a later step and can therefore be ignored
0128 
0129             %Recalculate the bounds for the new reaction. This can be
0130             %problematic since the scale of the bounds may change
0131             %dramatically. Let the most constraining reaction determine
0132             %the new bound
0133             lb1=reducedModel.lb(involvedRxns(1));
0134             lb2=reducedModel.lb(involvedRxns(2));
0135             ub1=reducedModel.ub(involvedRxns(1));
0136             ub2=reducedModel.ub(involvedRxns(2));
0137 
0138             if lb2~=-inf
0139                 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio);
0140             end
0141             if ub2~=inf
0142                 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio);
0143             end
0144             
0145             %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev
0146             reducedModel.rev(involvedRxns(1)) = reducedModel.rev(involvedRxns(1))*reducedModel.rev(involvedRxns(2));%this is a way to do an "and" operation with 0 and 1 numbers
0147 
0148             %Then recalculate the objective coefficient. The resulting
0149             %coefficient is the weighted sum of the previous
0150             reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio;
0151 
0152             %store which reactions that have been merged
0153             rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1))));
0154             rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2))));
0155             grpId = max(groupIds(rxnInd1),groupIds(rxnInd2));
0156             if grpId == 0
0157                grpId = nextGroupId;
0158                nextGroupId = nextGroupId + 1;
0159             end
0160 
0161             if groupIds(rxnInd1) ~= grpId
0162                if groupIds(rxnInd1) == 0
0163                    %not merged before, just set the group id
0164                    groupIds(rxnInd1) = grpId;
0165                else
0166                    %merged before - all rxns with the same group id should be changed
0167                    groupIds(groupIds == groupIds(rxnInd1)) = grpId;
0168                end
0169             end
0170             if groupIds(rxnInd2) ~= grpId
0171                if groupIds(rxnInd2) == 0
0172                    %not merged before, just set the group id
0173                    groupIds(rxnInd2) = grpId;
0174                else
0175                    %merged before - all rxns with the same group id should be changed
0176                    groupIds(groupIds == groupIds(rxnInd2)) = grpId;
0177                end
0178             end
0179 
0180             %Iterate again
0181             mergedSome=true;
0182         end
0183     end
0184 
0185     %All possible reactions merged
0186     if mergedSome==false
0187         break;
0188     end
0189 
0190     %Now delete all reactions that involve no metabolites
0191     I=find(sum(reducedModel.S~=0,1)==0);
0192 
0193     %Remove reactions
0194     reducedModel=removeReactions(reducedModel,I);
0195 
0196     %Remove metabolites
0197     notInUse=sum(reducedModel.S~=0,2)==0;
0198     reducedModel=removeMets(reducedModel,notInUse);
0199 end
0200     
0201 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1)
0202     model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1;
0203     %swap the bounds
0204     ub = model1.ub(rxnInd);
0205     model1.ub(rxnInd) = -model1.lb(rxnInd);
0206     model1.lb(rxnInd) = -ub;
0207     %flip the objective
0208     model1.c(rxnInd) = -model1.c(rxnInd);
0209     
0210     %now take care of the reversedRxns - if this is a group, reverse all of the
0211     %reactions in the group in the reversedRxns index - they will all be reversed at the
0212     %same time since they are the same rxn.
0213     rxnIndices = rxnInd;
0214     if groupIds1(rxnInd) > 0
0215         rxnIndices = find(groupIds1 == groupIds1(rxnInd));
0216     end
0217     reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices);
0218 end
0219 end