Home > RAVEN > sortModel.m

sortModel

PURPOSE ^

sortModel

SYNOPSIS ^

function model=sortModel(model,sortReversible,sortMetName,sortReactionOrder)

DESCRIPTION ^

 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
                     (opt, default true)
   sortMetName       sort the metabolite names in the equation, also uses
                     compartment abbreviation (opt, 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 (opt, default false)

   model             an updated model structure

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

   Rasmus Agren, 2013-02-06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005