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         if ~iscell(model.subSystems{i})
0069             model.subSystems{i} = {model.subSystems{i}};
0070         end    
0071         subsystemsConcatenated{i,1}=strjoin(model.subSystems{i,1},';');
0072         if ~isempty(model.subSystems{i,1})
0073             for j=1:numel(model.subSystems{i,1})
0074                 subsystemsUnique{numel(subsystemsUnique)+1,1}=model.subSystems{i,1}{1,j};
0075             end
0076         end
0077     end
0078     subsystemsUnique=unique(subsystemsUnique);
0079     for i=1:numel(subsystemsUnique)
0080         %Get all reactions for that subsystem
0081         rxns=find(~cellfun(@isempty,regexp(subsystemsConcatenated,subsystemsUnique(i))));
0082         
0083         %Temporarily ignore large subsystems because of inefficient
0084         %implementation
0085         if numel(rxns)<2 || numel(rxns)>250
0086             continue;
0087         end
0088         
0089         nRxns=numel(rxns);
0090         revRxns=rxns(model.rev(rxns)~=0);
0091         
0092         %This variable will hold the current reversibility directions of
0093         %the reversible reactions. 1 means the same direction as in the
0094         %original model and -1 means the opposite direction.
0095         oldRev=ones(numel(revRxns),1);
0096         
0097         %The problem could probably be solved analytically but a simple
0098         %random method is implemented here instead. Two reactions are
0099         %chosen randomly and their positions are switched. A score is
0100         %calculated based on the number of metabolites that are produced
0101         %before they are consumed. If the perturbed model has a better or
0102         %equal score than the original the reaction order is switched. If
0103         %no increase in score has been seen after 1000*rxnsInSubsystem then
0104         %the optimization is terminated
0105         
0106         rxnOrder=1:nRxns;
0107         oldScore=-inf;
0108         counter=0;
0109         firstIter=true;
0110         while 1==1
0111             counter=counter+1;
0112             if counter==100*nRxns
0113                 break;
0114             end
0115             
0116             newRxnOrder=rxnOrder;
0117             rev=oldRev;
0118             
0119             if firstIter==false
0120                 y=randperm(nRxns,2);
0121                 
0122                 %Switch the order
0123                 newRxnOrder(y(1))=rxnOrder(y(2));
0124                 newRxnOrder(y(2))=rxnOrder(y(1));
0125                 
0126                 %With a 50% chance, also switch the reversibility of one of
0127                 %the reactions
0128                 if rand()>0.5 && numel(rev)>1
0129                     n=randperm(numel(rev),1);
0130                     rev(n)=rev(n)*-1;
0131                 end
0132             end
0133             firstIter=false;
0134             
0135             tempS=model.S;
0136             
0137             %Switch the directionalities
0138             for j=1:numel(rev)
0139                 if rev(j)==-1
0140                     tempS(:,revRxns(j))=tempS(:,revRxns(j)).*-1;
0141                 end
0142             end
0143             
0144             %Get the metabolites that are involved and when they are
0145             %produced/consumed
0146             s=tempS(:,newRxnOrder);
0147             
0148             %Remove mets that aren't used in both directions
0149             s=s(any(s,2),:);
0150             
0151             %Add so that all mets are produced and consumed in the end
0152             s=[s ones(size(s,1),1) ones(size(s,1),1)*-1];
0153             
0154             %For each metabolite, find the reaction where it's first
0155             %produced and the reaction where it's first consumed
0156             s1=s>0;
0157             r1=arrayfun(@(x) find(s1(x,:),1,'first'),1:size(s1,1));
0158             s2=s<0;
0159             r2=arrayfun(@(x) find(s2(x,:),1,'first'),1:size(s2,1));
0160             
0161             score=sum(r1<r2);
0162             
0163             if score>=oldScore
0164                 if score>oldScore
0165                     counter=0;
0166                 end
0167                 oldScore=score;
0168                 oldRev=rev;
0169                 rxnOrder=newRxnOrder;
0170             end
0171         end
0172         
0173         %Update the model for this subsystem
0174         for j=1:numel(oldRev)
0175             if oldRev(j)==-1
0176                 model.S(:,revRxns(j))=model.S(:,revRxns(j)).*-1;
0177             end
0178         end
0179         order=1:numel(model.rxns);
0180         order(rxns)=rxns(rxnOrder);
0181         model=permuteModel(model, order, 'rxns');
0182     end
0183 end
0184 end