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         pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns));
0066         neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns));
0067         
0068                 
0069         %Check so that one or both of the reactions haven't been merged
0070         %already
0071         if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1
0072             %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible,
0073             %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise
0074             if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ...
0075                (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0)
0076                 involvedRxns = flip(involvedRxns);
0077             end
0078             
0079             %first make sure the first reaction is producing the metabolite
0080             if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0
0081                 %it is not producing the metabolite - fix that
0082                 %first choice: use the second reaction as producer if it is producing
0083                 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0
0084                     involvedRxns = flip(involvedRxns);%make the second the first
0085                 else
0086                     %now we know that the second reaction is not producing, so we can safely try to make the first a producer
0087                     if reducedModel.rev(involvedRxns(1)) == 1
0088                         [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns);
0089                     else %ok, finally try to flip the second reaction
0090                         if reducedModel.rev(involvedRxns(2)) == 1
0091                             [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns);
0092                             involvedRxns = flip(involvedRxns);%make the second the first
0093                         else
0094                             error('We should never end up here!');
0095                         end
0096                     end
0097                 end
0098             end
0099             %Now, make sure the second rxn is a consumer
0100             if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0
0101                 if reducedModel.rev(involvedRxns(2)) == 1
0102                     [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns);
0103                 else
0104                     error('We should never end up here!');
0105                 end
0106             end
0107 
0108             %Calculate how many times the second reaction has to be
0109             %multiplied before being merged with the first
0110             stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2)));
0111 
0112             %Add the second to the first
0113             reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio;
0114 
0115             %Clear the second reaction
0116             reducedModel.S(:,involvedRxns(2))=0;
0117 
0118             %This is to prevent numerical issues. It should be 0
0119             %already
0120             reducedModel.S(twoNonZero(i),involvedRxns(1))=0;
0121 
0122             %At this point the second reaction is certain to be deleted
0123             %in a later step and can therefore be ignored
0124 
0125             %Recalculate the bounds for the new reaction. This can be
0126             %problematic since the scale of the bounds may change
0127             %dramatically. Let the most constraining reaction determine
0128             %the new bound
0129             lb1=reducedModel.lb(involvedRxns(1));
0130             lb2=reducedModel.lb(involvedRxns(2));
0131             ub1=reducedModel.ub(involvedRxns(1));
0132             ub2=reducedModel.ub(involvedRxns(2));
0133 
0134             if lb2~=-inf
0135                 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio);
0136             end
0137             if ub2~=inf
0138                 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio);
0139             end
0140             
0141             %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev
0142             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
0143 
0144             %Then recalculate the objective coefficient. The resulting
0145             %coefficient is the weighted sum of the previous
0146             reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio;
0147 
0148             %store which reactions that have been merged
0149             rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1))));
0150             rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2))));
0151             grpId = max(groupIds(rxnInd1),groupIds(rxnInd2));
0152             if grpId == 0
0153                grpId = nextGroupId;
0154                nextGroupId = nextGroupId + 1;
0155             end
0156 
0157             if groupIds(rxnInd1) ~= grpId
0158                if groupIds(rxnInd1) == 0
0159                    %not merged before, just set the group id
0160                    groupIds(rxnInd1) = grpId;
0161                else
0162                    %merged before - all rxns with the same group id should be changed
0163                    groupIds(groupIds == groupIds(rxnInd1)) = grpId;
0164                end
0165             end
0166             if groupIds(rxnInd2) ~= grpId
0167                if groupIds(rxnInd2) == 0
0168                    %not merged before, just set the group id
0169                    groupIds(rxnInd2) = grpId;
0170                else
0171                    %merged before - all rxns with the same group id should be changed
0172                    groupIds(groupIds == groupIds(rxnInd2)) = grpId;
0173                end
0174             end
0175 
0176             %Iterate again
0177             mergedSome=true;
0178         end
0179     end
0180 
0181     %All possible reactions merged
0182     if mergedSome==false
0183         break;
0184     end
0185 
0186     %Now delete all reactions that involve no metabolites
0187     I=find(sum(reducedModel.S~=0,1)==0);
0188 
0189     %Remove reactions
0190     reducedModel=removeReactions(reducedModel,I);
0191 
0192     %Remove metabolites
0193     notInUse=sum(reducedModel.S~=0,2)==0;
0194     reducedModel=removeMets(reducedModel,notInUse);
0195 end
0196     
0197 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1)
0198     model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1;
0199     %swap the bounds
0200     ub = model1.ub(rxnInd);
0201     model1.ub(rxnInd) = -model1.lb(rxnInd);
0202     model1.lb(rxnInd) = -ub;
0203     %flip the objective
0204     model1.c(rxnInd) = -model1.c(rxnInd);
0205     
0206     %now take care of the reversedRxns - if this is a group, reverse all of the
0207     %reactions in the group in the reversedRxns index - they will all be reversed at the
0208     %same time since they are the same rxn.
0209     rxnIndices = rxnInd;
0210     if groupIds1(rxnInd) > 0
0211         rxnIndices = find(groupIds1 == groupIds1(rxnInd));
0212     end
0213     reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices);
0214 end
0215 end