expandModel

 expandModel
   Expands a model which uses several gene associations for one reaction.
   Each such reaction is split into several reactions, each under the control
   of only one gene.
  
 Input:
   model       model structure
 
 Output:
   newModel    model structure with separate reactions for iso-enzymes, where
               the reaction ids are renamed as to id_EXP_1, id_EXP_2, etc. 
   rxnToCheck  cell array with original reaction identifiers for those
               that contained nested and/or-relationships in grRules.

   NOTE: grRules strings that involve nested expressions of 'and' and 'or'
         might not be parsed correctly if they are not standardized (if the
         standardizeGrRules functions was not first run on the model). For
         those reactions, it is therefore advisable to inspect the reactions in
         rxnToCheck to confirm correct model expansion.

 Usage: [newModel, rxnToCheck]=expandModel(model)

0001 function [newModel, rxnToCheck]=expandModel(model)
0002 % expandModel
0003 %   Expands a model which uses several gene associations for one reaction.
0004 %   Each such reaction is split into several reactions, each under the control
0005 %   of only one gene.
0006 %
0007 % Input:
0008 %   model       model structure
0009 %
0010 % Output:
0011 %   newModel    model structure with separate reactions for iso-enzymes, where
0012 %               the reaction ids are renamed as to id_EXP_1, id_EXP_2, etc.
0013 %   rxnToCheck  cell array with original reaction identifiers for those
0014 %               that contained nested and/or-relationships in grRules.
0015 %
0016 %   NOTE: grRules strings that involve nested expressions of 'and' and 'or'
0017 %         might not be parsed correctly if they are not standardized (if the
0018 %         standardizeGrRules functions was not first run on the model). For
0019 %         those reactions, it is therefore advisable to inspect the reactions in
0020 %         rxnToCheck to confirm correct model expansion.
0021 %
0022 % Usage: [newModel, rxnToCheck]=expandModel(model)
0023 
0024 %Check how many reactions we will create (the number of or:s in the GPRs).
0025 %This way, we can preallocate all fields and save much computation time
0026 
0027 numOrs = count(model.grRules, ' or ');
0028 toAdd = sum(numOrs);
0029 prevNumRxns = length(model.rxns);
0030 rxnToCheck={};
0031 if toAdd > 0
0032     %Calculate indices to expand
0033     %For example, if a reaction with index x has 2 or:s, meaning it has 3
0034     %reactions after the split, we should add two copies of this reaction
0035     %For fields that should just be copied to the new reactions, we just keep
0036     %track of that there are two copies, i.e., we add x x to this vector.
0037     %That is exactly what repelem does for us.
0038     cpyIndices = repelem(1:prevNumRxns, numOrs);
0039     
0040     %Copy all fields that should just be copied
0041     model.S=[model.S model.S(:,cpyIndices)];
0042     model.rxnNames=[model.rxnNames;model.rxnNames(cpyIndices)];
0043     model.lb=[model.lb;model.lb(cpyIndices)];
0044     model.ub=[model.ub;model.ub(cpyIndices)];
0045     model.rev=[model.rev;model.rev(cpyIndices)];
0046     model.c=[model.c;model.c(cpyIndices)];
0047     if isfield(model,'subSystems')
0048         model.subSystems=[model.subSystems;model.subSystems(cpyIndices)];
0049     end
0050     if isfield(model,'eccodes')
0051         model.eccodes=[model.eccodes;model.eccodes(cpyIndices)];
0052     end
0053     if isfield(model,'equations')
0054         model.equations=[model.equations;model.equations(cpyIndices)];
0055     end
0056     if isfield(model,'rxnMiriams')
0057         model.rxnMiriams=[model.rxnMiriams;model.rxnMiriams(cpyIndices)];
0058     end
0059     if isfield(model,'rxnComps')
0060         model.rxnComps=[model.rxnComps;model.rxnComps(cpyIndices)];
0061     end
0062     if isfield(model,'rxnFrom')
0063         model.rxnFrom=[model.rxnFrom;model.rxnFrom(cpyIndices)];
0064     end
0065     if isfield(model,'rxnNotes')
0066         model.rxnNotes=[model.rxnNotes;model.rxnNotes(cpyIndices)];
0067     end
0068     if isfield(model,'rxnReferences')
0069         model.rxnReferences=[model.rxnReferences;model.rxnReferences(cpyIndices)];
0070     end
0071     if isfield(model,'rxnConfidenceScores')
0072         model.rxnConfidenceScores=[model.rxnConfidenceScores;model.rxnConfidenceScores(cpyIndices)];
0073     end
0074     if isfield(model,'rxnDeltaG')
0075         model.rxnDeltaG=[model.rxnDeltaG;model.rxnDeltaG(cpyIndices)];
0076     end
0077     
0078     %now expand the more complex fields - will be filled in later
0079     model.rxns=[model.rxns;cell(toAdd,1)];
0080     model.grRules=[model.grRules;cell(toAdd,1)];
0081     model.rxnGeneMat=[model.rxnGeneMat;sparse(toAdd,size(model.rxnGeneMat,2))];
0082     
0083     %Loop throught those reactions and fill in the expanded data
0084     nextIndex = prevNumRxns + 1;
0085     for i=1:prevNumRxns
0086         if (numOrs(i) > 0)
0087             %Check that it doesn't contain nested 'and' and 'or' relations and
0088             %print a warning if it does
0089             if ~isempty(strfind(model.grRules{i},' and '))
0090                 rxnToCheck{end+1,1}=model.rxns{i};
0091             end
0092 
0093             %Get rid of all '(' and ')' since I'm not looking at complex stuff
0094             %anyways
0095             geneString=model.grRules{i};
0096             geneString=strrep(geneString,'(','');
0097             geneString=strrep(geneString,')','');
0098             geneString=strrep(geneString,' or ',';');
0099 
0100             %Split the string into gene names
0101             geneNames=regexp(geneString,';','split');
0102 
0103             %Update the reaction to only use the first gene
0104             model.grRules{i}=['(' geneNames{1} ')'];
0105             %Find the gene in the gene list If ' and ' relationship, first
0106             %split the genes
0107             model.rxnGeneMat(i,:)=0;
0108             if ~isempty(strfind(geneNames(1),' and '))
0109                 andGenes=regexp(geneNames{1},' and ','split');
0110                 model.rxnGeneMat(i,ismember(model.genes,andGenes)) = 1;
0111             else
0112                 [~, index]=ismember(geneNames(1),model.genes);
0113                 model.rxnGeneMat(i,index)=1;
0114             end
0115 
0116             %Insert the reactions at the end of the model and without
0117             %allocating space. This is not nice, but ok for now
0118             for j=2:numel(geneNames)
0119                 ind = nextIndex+j-2;
0120                 model.rxns{ind}=[model.rxns{i} '_EXP_' num2str(j)];
0121                 
0122                 model.grRules{ind}=['(' geneNames{j} ')'];
0123 
0124                 if ~isempty(strfind(geneNames(j),' and '))
0125                     andGenes=regexp(geneNames{j},' and ','split');
0126                     model.rxnGeneMat(ind,ismember(model.genes,andGenes)) = 1;
0127                 else
0128                     model.rxnGeneMat(ind,ismember(model.genes,geneNames(j))) = 1;
0129                 end
0130             end
0131             model.rxns{i}=[model.rxns{i}, '_EXP_1'];
0132             nextIndex = nextIndex + numOrs(i);
0133         end        
0134     end
0135     newModel=model;
0136 else
0137     %There are no reactions to expand, return the model as is
0138     newModel=model;
0139 end
0140 
0141 %Fix grRules and reconstruct rxnGeneMat
0142 [grRules,rxnGeneMat] = standardizeGrRules(newModel,true);
0143 newModel.grRules = grRules;
0144 newModel.rxnGeneMat = rxnGeneMat;
0145 end