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