sortModel

 sortModel
   Sorts a model based on metabolite names and compartments

   model             a model structure
   sortReversible    sorts the reversible reactions so the the metabolite
                     that is first in lexiographical order is a reactant
                     (optional, default true)
   sortMetName       sort the metabolite names in the equation, also uses
                     compartment abbreviation (optional, default false)
   sortReactionOrder sorts the reaction order within each subsystem so that
                     reactions consuming some metabolite comes efter
                     reactions producing it. This overrides the
                     sortReversible option and reactions are sorted so that
                     the production direction matches the consumption
                     direction (optional, default false)

   model             an updated model structure

 Usage: model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)

0001 function model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)
0002 % sortModel
0003 %   Sorts a model based on metabolite names and compartments
0004 %
0005 %   model             a model structure
0006 %   sortReversible    sorts the reversible reactions so the the metabolite
0007 %                     that is first in lexiographical order is a reactant
0008 %                     (optional, default true)
0009 %   sortMetName       sort the metabolite names in the equation, also uses
0010 %                     compartment abbreviation (optional, default false)
0011 %   sortReactionOrder sorts the reaction order within each subsystem so that
0012 %                     reactions consuming some metabolite comes efter
0013 %                     reactions producing it. This overrides the
0014 %                     sortReversible option and reactions are sorted so that
0015 %                     the production direction matches the consumption
0016 %                     direction (optional, default false)
0017 %
0018 %   model             an updated model structure
0019 %
0020 % Usage: model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)
0021 
0022 if nargin<2
0023     sortReversible=true;
0024 end
0025 if nargin<3
0026     sortMetName=false;
0027 end
0028 if nargin<4
0029     sortReactionOrder=false;
0030 end
0031 
0032 if sortMetName==true
0033     %Assuming that metComps are the indexes. Should be changed at one point
0034     [~, metIndexes]=sort(strcat(model.metNames,'[',model.comps(model.metComps),']'));
0035     model=permuteModel(model,metIndexes,'mets');
0036 end
0037 
0038 if sortReversible==true && sortReactionOrder==false
0039     %Get all reversible reactions
0040     revIndexes=find(model.rev);
0041     
0042     %Loop through them
0043     for i=1:numel(revIndexes)
0044         %Create a cell array with all the metabolite names
0045         mets=find(model.S(:,revIndexes(i)));
0046         metNames=strcat(model.metNames(mets),model.comps(model.metComps(mets)));
0047         
0048         if iscellstr(metNames)
0049             [~, indexes]=sort(metNames);
0050             
0051             if model.S(mets(indexes(1)),revIndexes(i))>0
0052                 model.S(:,revIndexes(i))=model.S(:,revIndexes(i))*-1;
0053             end
0054         end
0055     end
0056 end
0057 
0058 if sortReactionOrder==true
0059     %Check if the model has sub-systems, otherwise throw an error
0060     if ~isfield(model,'subSystems')
0061         EM='The model must contain a subSystems field in order to sort reaction order';
0062         dispEM(EM);
0063     end
0064     
0065     subsystemsUnique='';
0066     subsystemsConcatenated='';
0067     for i=1:numel(model.subSystems)
0068         subsystemsConcatenated{i,1}=strjoin(model.subSystems{i,1},';');
0069         if ~isempty(model.subSystems{i,1})
0070             for j=1:numel(model.subSystems{i,1})
0071                 subsystemsUnique{numel(subsystemsUnique)+1,1}=model.subSystems{i,1}{1,j};
0072             end
0073         end
0074     end
0075     subsystemsUnique=unique(subsystemsUnique);
0076     for i=1:numel(subsystemsUnique)
0077         %Get all reactions for that subsystem
0078         rxns=find(~cellfun(@isempty,regexp(subsystemsConcatenated,subsystemsUnique(i))));
0079         
0080         %Temporarily ignore large subsystems because of inefficient
0081         %implementation
0082         if numel(rxns)<2 || numel(rxns)>250
0083             continue;
0084         end
0085         
0086         nRxns=numel(rxns);
0087         revRxns=rxns(model.rev(rxns)~=0);
0088         
0089         %This variable will hold the current reversibility directions of
0090         %the reversible reactions. 1 means the same direction as in the
0091         %original model and -1 means the opposite direction.
0092         oldRev=ones(numel(revRxns),1);
0093         
0094         %The problem could probably be solved analytically but a simple
0095         %random method is implemented here instead. Two reactions are
0096         %chosen randomly and their positions are switched. A score is
0097         %calculated based on the number of metabolites that are produced
0098         %before they are consumed. If the perturbed model has a better or
0099         %equal score than the original the reaction order is switched. If
0100         %no increase in score has been seen after 1000*rxnsInSubsystem then
0101         %the optimization is terminated
0102         
0103         rxnOrder=1:nRxns;
0104         oldScore=-inf;
0105         counter=0;
0106         firstIter=true;
0107         while 1==1
0108             counter=counter+1;
0109             if counter==100*nRxns
0110                 break;
0111             end
0112             
0113             newRxnOrder=rxnOrder;
0114             rev=oldRev;
0115             
0116             if firstIter==false
0117                 y=randperm(nRxns,2);
0118                 
0119                 %Switch the order
0120                 newRxnOrder(y(1))=rxnOrder(y(2));
0121                 newRxnOrder(y(2))=rxnOrder(y(1));
0122                 
0123                 %With a 50% chance, also switch the reversibility of one of
0124                 %the reactions
0125                 if rand()>0.5 && numel(rev)>1
0126                     n=randperm(numel(rev),1);
0127                     rev(n)=rev(n)*-1;
0128                 end
0129             end
0130             firstIter=false;
0131             
0132             tempS=model.S;
0133             
0134             %Switch the directionalities
0135             for j=1:numel(rev)
0136                 if rev(j)==-1
0137                     tempS(:,revRxns(j))=tempS(:,revRxns(j)).*-1;
0138                 end
0139             end
0140             
0141             %Get the metabolites that are involved and when they are
0142             %produced/consumed
0143             s=tempS(:,newRxnOrder);
0144             
0145             %Remove mets that aren't used in both directions
0146             s=s(any(s,2),:);
0147             
0148             %Add so that all mets are produced and consumed in the end
0149             s=[s ones(size(s,1),1) ones(size(s,1),1)*-1];
0150             
0151             %For each metabolite, find the reaction where it's first
0152             %produced and the reaction where it's first consumed
0153             s1=s>0;
0154             r1=arrayfun(@(x) find(s1(x,:),1,'first'),1:size(s1,1));
0155             s2=s<0;
0156             r2=arrayfun(@(x) find(s2(x,:),1,'first'),1:size(s2,1));
0157             
0158             score=sum(r1<r2);
0159             
0160             if score>=oldScore
0161                 if score>oldScore
0162                     counter=0;
0163                 end
0164                 oldScore=score;
0165                 oldRev=rev;
0166                 rxnOrder=newRxnOrder;
0167             end
0168         end
0169         
0170         %Update the model for this subsystem
0171         for j=1:numel(oldRev)
0172             if oldRev(j)==-1
0173                 model.S(:,revRxns(j))=model.S(:,revRxns(j)).*-1;
0174             end
0175         end
0176         order=1:numel(model.rxns);
0177         order(rxns)=rxns(rxnOrder);
0178         model=permuteModel(model, order, 'rxns');
0179     end
0180 end
0181 end