mergeModels

 mergeModels
   Merges models into one model structure. Reactions are added without any
   checks, so duplicate reactions might appear. Metabolites are matched by
   their name and compartment (metaboliteName[comp]), while genes are
   matched by their name.

   models          a cell array with model structures
   metParam        string specifying whether to refer to metabolite name
                   (metNames) or ID (mets) for matching (default, metNames)
   supressWarnings true if warnings should be supressed (optional, default
                   false)

   model     a model structure with the merged model. Follows the structure
             of normal models but also has 'rxnFrom/metFrom/geneFrom' fields
             to indicate from which model each reaction/metabolite/gene was
             taken

 Usage: model=mergeModels(models)

0001 function model=mergeModels(models,metParam,supressWarnings)
0002 % mergeModels
0003 %   Merges models into one model structure. Reactions are added without any
0004 %   checks, so duplicate reactions might appear. Metabolites are matched by
0005 %   their name and compartment (metaboliteName[comp]), while genes are
0006 %   matched by their name.
0007 %
0008 %   models          a cell array with model structures
0009 %   metParam        string specifying whether to refer to metabolite name
0010 %                   (metNames) or ID (mets) for matching (default, metNames)
0011 %   supressWarnings true if warnings should be supressed (optional, default
0012 %                   false)
0013 %
0014 %   model     a model structure with the merged model. Follows the structure
0015 %             of normal models but also has 'rxnFrom/metFrom/geneFrom' fields
0016 %             to indicate from which model each reaction/metabolite/gene was
0017 %             taken
0018 %
0019 % Usage: model=mergeModels(models)
0020 
0021 %Just return the model
0022 if numel(models)<=1
0023     model=models{1};
0024     return;
0025 end
0026 
0027 if nargin<2
0028     metParam='metNames';
0029 else
0030     metParam=char(metParam);
0031 end
0032 
0033 if nargin<3
0034     supressWarnings=false;
0035 end
0036 
0037 %Add new functionality in the order specified in models
0038 model=models{1};
0039 model.id='MERGED';
0040 model.name='';
0041 
0042 model.rxnFrom=cell(numel(models{1}.rxns),1);
0043 model.rxnFrom(:)={models{1}.id};
0044 model.metFrom=cell(numel(models{1}.mets),1);
0045 model.metFrom(:)={models{1}.id};
0046 if isfield(models{1},'genes')
0047     model.geneFrom=cell(numel(models{1}.genes),1);
0048     model.geneFrom(:)={models{1}.id};
0049 end
0050 
0051 if isfield(model,'equations')
0052     model=rmfield(model,'equations');
0053 end
0054 
0055 for i=2:numel(models)
0056     %Add the model id to the rxn id id it already exists in the model (id
0057     %have to be unique) This is because it makes a '[]' string if no new
0058     %reactions
0059     if ~isempty(models{i}.rxns)
0060         I=ismember(models{i}.rxns,model.rxns);
0061         models{i}.rxns(I)=strcat(models{i}.rxns(I),['_' models{i}.id]);
0062     end
0063     
0064     %Make sure that there are no conflicting reaction ids
0065     [~, ~, conflicting]=intersect(model.rxns,models{i}.rxns);
0066     
0067     if ~isempty(conflicting)
0068         printString=cell(numel(conflicting),1);
0069         for j=1:numel(conflicting)
0070             printString{j}=['Old: ' models{i}.rxns{conflicting(j)} ' New: ' models{i}.rxns{conflicting(j)} '_' models{i}.id];
0071             models{i}.rxns{conflicting(j)}=[models{i}.rxns{conflicting(j)} '_' models{i}.id];
0072         end
0073         if supressWarnings==false
0074             EM=['The following reaction IDs in ' models{i}.id ' are already present in the model and were renamed:'];
0075             dispEM(EM,false,printString);
0076             fprintf('\n');
0077         end
0078     end
0079     
0080     %Add all static stuff
0081     rxnFrom=cell(numel(models{i}.rxns),1);
0082     rxnFrom(:)={models{i}.id};
0083     model.rxnFrom=[model.rxnFrom;rxnFrom];
0084     model.rxns=[model.rxns;models{i}.rxns];
0085     model.rxnNames=[model.rxnNames;models{i}.rxnNames];
0086     model.lb=[model.lb;models{i}.lb];
0087     model.ub=[model.ub;models{i}.ub];
0088     model.c=[model.c;models{i}.c];
0089     model.rev=[model.rev;models{i}.rev];
0090     
0091     if isfield(models{i},'subSystems')
0092         if isfield(model,'subSystems')
0093             model.subSystems=[model.subSystems;models{i}.subSystems];
0094         else
0095             emptySubSystem=cell(numel(model.rxns)-numel(models{i}.rxns),1);
0096             emptySubSystem(:)={''};
0097             emptySubSystem=cellfun(@(x) cell(0,0),emptySubSystem,'UniformOutput',false);
0098             model.subSystems=[emptySubSystem;models{i}.subSystems];
0099         end
0100     else
0101         if isfield(model,'subSystems')
0102             emptySubSystem=cell(numel(models{i}.rxns),1);
0103             emptySubSystem(:)={''};
0104             emptySubSystem=cellfun(@(x) cell(0,0),emptySubSystem,'UniformOutput',false);
0105             model.subSystems=[model.subSystems;emptySubSystem];
0106         end
0107     end
0108     
0109     if isfield(models{i},'eccodes')
0110         if isfield(model,'eccodes')
0111             model.eccodes=[model.eccodes;models{i}.eccodes];
0112         else
0113             emptyEC=cell(numel(model.rxns)-numel(models{i}.rxns),1);
0114             emptyEC(:)={''};
0115             model.eccodes=[emptyEC;models{i}.eccodes];
0116         end
0117     else
0118         if isfield(model,'eccodes')
0119             emptyEC=cell(numel(models{i}.rxns),1);
0120             emptyEC(:)={''};
0121             model.eccodes=[model.eccodes;emptyEC];
0122         end
0123     end
0124     
0125     if isfield(models{i},'rxnMiriams')
0126         if isfield(model,'rxnMiriams')
0127             model.rxnMiriams=[model.rxnMiriams;models{i}.rxnMiriams];
0128         else
0129             model.rxnMiriams=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnMiriams];
0130         end
0131     else
0132         if isfield(model,'rxnMiriams')
0133             model.rxnMiriams=[model.rxnMiriams;cell(numel(models{i}.rxns),1)];
0134         end
0135     end
0136     
0137     if isfield(models{i},'rxnNotes')
0138         if isfield(model,'rxnNotes')
0139             model.rxnNotes=[model.rxnNotes;models{i}.rxnNotes];
0140         else
0141             emptyNotes=cell(numel(model.rxns)-numel(models{i}.rxns),1);
0142             emptyNotes(:)={''};
0143             model.rxnNotes=[emptyNotes;models{i}.rxnNotes];
0144         end
0145     else
0146         if isfield(model,'rxnNotes')
0147             emptyNotes=cell(numel(models{i}.rxns),1);
0148             emptyNotes(:)={''};
0149             model.rxnNotes=[model.rxnNotes;emptyNotes];
0150         end
0151     end
0152     
0153     if isfield(models{i},'rxnReferences')
0154         if isfield(model,'rxnReferences')
0155             model.rxnReferences=[model.rxnReferences;models{i}.rxnReferences];
0156         else
0157             emptyReferences=cell(numel(model.rxns)-numel(models{i}.rxns),1);
0158             emptyReferences(:)={''};
0159             model.rxnReferences=[emptyReferences;models{i}.rxnReferences];
0160         end
0161     else
0162         if isfield(model,'rxnReferences')
0163             emptyReferences=cell(numel(models{i}.rxns),1);
0164             emptyReferences(:)={''};
0165             model.rxnReferences=[model.rxnReferences;emptyReferences];
0166         end
0167     end
0168     
0169     if isfield(models{i},'rxnConfidenceScores')
0170         if isfield(model,'rxnConfidenceScores')
0171             model.rxnConfidenceScores=[model.rxnConfidenceScores;models{i}.rxnConfidenceScores];
0172         else
0173             model.rxnConfidenceScores=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnConfidenceScores];
0174         end
0175     else
0176         if isfield(model,'rxnConfidenceScores')
0177             model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(models{i}.rxns),1)];
0178         end
0179     end
0180 
0181     if isfield(models{i},'rxnDeltaG')
0182         if isfield(model,'rxnDeltaG')
0183             model.rxnDeltaG=[model.rxnDeltaG;models{i}.rxnDeltaG];
0184         else
0185             model.rxnDeltaG=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnDeltaG];
0186         end
0187     else
0188         if isfield(model,'rxnDeltaG')
0189             model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(models{i}.rxns),1)];
0190         end
0191     end
0192     
0193     if isfield(models{i},'rxnComps')
0194         if isfield(model,'rxnComps')
0195             model.rxnComps=[model.rxnComps;models{i}.rxnComps];
0196         else
0197             model.rxnComps=[ones(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnComps];
0198             fprintf('NOTE: One of the models does not contain compartment information for its reactions. All reactions in that model has been assigned to the first compartment\n');
0199         end
0200     else
0201         if isfield(model,'rxnComps')
0202             model.rxnComps=[model.rxnComps;ones(numel(models{i}.rxns),1)];
0203             fprintf('NOTE: One of the models does not contain compartment information for its reactions. All reactions in that model has been assigned to the first compartment\n');
0204         end
0205     end
0206     
0207     if isfield(models{i},'rxnScores')
0208         if isfield(model,'rxnScores')
0209             model.rxnScores=[model.rxnScores;models{i}.rxnScores];
0210         else
0211             emptyRS=zeros(numel(model.rxns)-numel(models{i}.rxns),1);
0212             model.rxnScores=[emptyRS;models{i}.rxnScores];
0213         end
0214     else
0215         if isfield(model,'rxnScores')
0216             emptyRS=zeros(numel(models{i}.rxns),1);
0217             model.rxnScores=[model.rxnScores;emptyRS];
0218         end
0219     end
0220     
0221     if isfield(models{i},'pwys')
0222         if isfield(model,'pwys')
0223             model.pwys=[model.pwys;models{i}.pwys];
0224         else
0225             model.pwys=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.pwys];
0226         end
0227     else
0228         if isfield(model,'pwys')
0229             model.pwys=[model.pwys;cell(numel(models{i}.rxns),1)];
0230         end
0231     end
0232 
0233     if strcmpi(metParam,'metNames')
0234     %Get the new metabolites from matching the models. Metabolites are said
0235     %to be the same if they share name and compartment id. This means that
0236     %metabolite IDs are not taken into account.
0237         
0238         oldMetComps=model.comps(model.metComps);
0239         oldMets=strcat(model.metNames,'[',oldMetComps,']');
0240         %This is because it makes a '[]' string if no new metabolites
0241         if ~isempty(models{i}.metNames)
0242             newMetComps=models{i}.comps(models{i}.metComps);
0243             newMets=strcat(models{i}.metNames,'[',newMetComps,']');
0244         else
0245             newMets={};
0246             newMetComps={};
0247         end
0248         tf=ismember(newMets,oldMets);
0249         metsToAdd=find(~tf);
0250 
0251     end
0252 
0253     if strcmpi(metParam,'mets')
0254     %Get the new metabolites from matching the models. Metabolites are matched by metabolite ID (model.mets).
0255 
0256         oldMetComps=model.comps(model.metComps);
0257         oldMets=model.mets;
0258     
0259         if ~isempty(models{i}.mets)
0260             newMetComps=models{i}.comps(models{i}.metComps);
0261             newMets=models{i}.mets;
0262         else
0263             newMets={};
0264             newMetComps={};
0265         end
0266         tf=ismember(newMets,oldMets);
0267         metsToAdd=find(~tf);
0268 
0269     end
0270     
0271     %First add the new metabolites Make sure that there are no conflicting
0272     %metabolite ids
0273     conflicting=ismember(models{i}.mets(metsToAdd),model.mets);
0274     
0275     conflicting=find(conflicting);
0276     
0277     if ~isempty(conflicting)
0278         printString=cell(numel(conflicting),1);
0279         for j=1:numel(conflicting)
0280             printString{j}=['Old: ' models{i}.mets{metsToAdd(conflicting(j))} ' New: ' models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
0281             models{i}.mets{metsToAdd(conflicting(j))}=[models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
0282         end
0283         if supressWarnings==false
0284             EM=['The following metabolite IDs in ' models{i}.id ' are already present in the model and were renamed:'];
0285             dispEM(EM,false,printString);
0286         end
0287     end
0288     
0289     %Add static info on the metabolites
0290     metFrom=cell(numel(metsToAdd),1);
0291     metFrom(:)={models{i}.id};
0292     model.metFrom=[model.metFrom;metFrom];
0293     model.mets=[model.mets;models{i}.mets(metsToAdd)];
0294     model.metNames=[model.metNames;models{i}.metNames(metsToAdd)];
0295     model.b=[model.b;zeros(numel(metsToAdd),size(model.b,2))];
0296     
0297     if isfield(model,'unconstrained')
0298         if isfield(models{i},'unconstrained')
0299             model.unconstrained=[model.unconstrained;models{i}.unconstrained(metsToAdd)];
0300         else
0301             model.unconstrained=[model.unconstrained;zeros(numel(metsToAdd),1)];
0302         end
0303     else
0304         if isfield(models{i},'unconstrained')
0305             model.unconstrained=[zeros(numel(model.mets),1);models{i}.unconstrained(metsToAdd)];
0306         end
0307     end
0308     
0309     %Only add extra info on new metabolites since it's a little tricky to
0310     %chose what to keep otherwise. Should change in the future
0311 
0312     if ~isempty(metsToAdd)
0313         if isfield(models{i},'inchis')
0314             if isfield(model,'inchis')
0315                 model.inchis=[model.inchis;models{i}.inchis(metsToAdd)];
0316             else
0317                 emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
0318                 emptyInchi(:)={''};
0319                 model.inchis=[emptyInchi;models{i}.inchis(metsToAdd)];
0320             end
0321         else
0322             if isfield(model,'inchis')
0323                 emptyInchi=cell(numel(metsToAdd),1);
0324                 emptyInchi(:)={''};
0325                 model.inchis=[model.inchis;emptyInchi];
0326             end
0327         end
0328 
0329         if isfield(models{i},'metSmiles')
0330             if isfield(model,'metSmiles')
0331                 model.metSmiles=[model.metSmiles;models{i}.metSmiles(metsToAdd)];
0332             else
0333                 emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
0334                 emptyInchi(:)={''};
0335                 model.metSmiles=[emptyInchi;models{i}.metSmiles(metsToAdd)];
0336             end
0337         else
0338             if isfield(model,'metSmiles')
0339                 emptyInchi=cell(numel(metsToAdd),1);
0340                 emptyInchi(:)={''};
0341                 model.metSmiles=[model.metSmiles;emptyInchi];
0342             end
0343         end
0344         
0345         if isfield(models{i},'metFormulas')
0346             if isfield(model,'metFormulas')
0347                 model.metFormulas=[model.metFormulas;models{i}.metFormulas(metsToAdd)];
0348             else
0349                 emptyMetFormulas=cell(numel(model.mets)-numel(metsToAdd),1);
0350                 emptyMetFormulas(:)={''};
0351                 model.metFormulas=[emptyMetFormulas;models{i}.metFormulas(metsToAdd)];
0352             end
0353         else
0354             if isfield(model,'metFormulas')
0355                 emptyMetFormulas=cell(numel(metsToAdd),1);
0356                 emptyMetFormulas(:)={''};
0357                 model.metFormulas=[model.metFormulas;emptyMetFormulas];
0358             end
0359         end
0360         
0361         if isfield(models{i},'metCharges')
0362             if isfield(model,'metCharges')
0363                 model.metCharges=[model.metCharges;models{i}.metCharges(metsToAdd)];
0364             else
0365                 emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
0366                 model.metCharges=[emptyMetCharge;models{i}.metCharges(metsToAdd)];
0367             end
0368         else
0369             if isfield(model,'metCharges')
0370                 emptyMetCharge=nan(numel(metsToAdd),1);
0371                 model.metCharges=[model.metCharges;emptyMetCharge];
0372             end
0373         end
0374 
0375         if isfield(models{i},'metDeltaG')
0376             if isfield(model,'metDeltaG')
0377                 model.metDeltaG=[model.metDeltaG;models{i}.metDeltaG(metsToAdd)];
0378             else
0379                 emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
0380                 model.metDeltaG=[emptyMetCharge;models{i}.metDeltaG(metsToAdd)];
0381             end
0382         else
0383             if isfield(model,'metDeltaG')
0384                 emptyMetCharge=nan(numel(metsToAdd),1);
0385                 model.metDeltaG=[model.metDeltaG;emptyMetCharge];
0386             end
0387         end
0388         
0389         if isfield(models{i},'metMiriams')
0390             if isfield(model,'metMiriams')
0391                 model.metMiriams=[model.metMiriams;models{i}.metMiriams(metsToAdd)];
0392             else
0393                 emptyMetMiriam=cell(numel(model.mets)-numel(metsToAdd),1);
0394                 model.metMiriams=[emptyMetMiriam;models{i}.metMiriams(metsToAdd)];
0395             end
0396         else
0397             if isfield(model,'metMiriams')
0398                 emptyMetMiriam=cell(numel(metsToAdd),1);
0399                 model.metMiriams=[model.metMiriams;emptyMetMiriam];
0400             end
0401         end
0402     end
0403     
0404     %Add if there are any new compartments and add those. This can change
0405     %the order of compartments and the corresponding indexes in
0406     %model.metComps.
0407     
0408     %Find overlapping and new compartments
0409     [overlap, oldIDs]=ismember(models{i}.comps,model.comps);
0410     overlap=find(overlap);
0411     
0412     %Add the new compartments if any
0413     if numel(overlap)~=numel(models{i}.compNames)
0414         compIndexes=oldIDs==0;
0415         
0416         %Make sure that there are no conflicting compartment ids
0417         [~, conflicting]=ismember(models{i}.compNames(compIndexes),model.compNames);
0418         if any(conflicting)
0419             EM=['The following compartment IDs in ' models{i}.id ' are already present in the model but with another name. They have to be renamed'];
0420             dispEM(EM,true,model.comps(conflicting));
0421         end
0422         
0423         %It's ok to add duplicate name, but not duplicate IDs
0424         model.compNames=[model.compNames; models{i}.compNames(compIndexes)];
0425         model.comps=[model.comps; models{i}.comps(compIndexes)];
0426         if isfield(model,'compOutside')
0427             if isfield(models{i},'compOutside')
0428                 model.compOutside=[model.compOutside; models{i}.compOutside(compIndexes)];
0429             else
0430                 %This is if not all models have the field
0431                 padding=cell(sum(compIndexes),1);
0432                 padding(:)={''};
0433                 model.compOutside=[model.compOutside;padding];
0434             end
0435         end
0436         if isfield(model,'compMiriams')
0437             if isfield(models{i},'compMiriams')
0438                 model.compMiriams=[model.compMiriams; models{i}.compMiriams(compIndexes)];
0439             else
0440                 %This is if not all models have the field
0441                 model.compMiriams=[model.compMiriams;cell(sum(compIndexes),1)];
0442             end
0443         end
0444     end
0445     
0446     %Only add new comp info on the un-matched metabolites since the old
0447     %ones will be mapped to the existing list anyways
0448     [I, J]=ismember(newMetComps(metsToAdd),model.comps);
0449     %Just a check
0450     if ~all(I)
0451         EM='There was an unexpected error in matching compartments';
0452         dispEM(EM);
0453     end
0454     model.metComps=[model.metComps;J];
0455      
0456     %Create the new stoichiometric matrix
0457     model.S=[model.S;sparse(numel(metsToAdd),size(model.S,2))];
0458     
0459 
0460     if strcmpi(metParam,'metNames')
0461         %Rematch metabolite names. Not the most clever way to do it maybe
0462         allMets=strcat(model.metNames,'[',model.comps(model.metComps),']');
0463         [~, J]=ismember(newMets,allMets);
0464     end
0465 
0466     if strcmpi(metParam,'mets')
0467         %Rematch metabolite by IDs and add unique new metabolites
0468         allMets=model.mets;
0469         uniqueNewMets = setdiff(newMets,oldMets);
0470         allMets(end+1:end+numel(uniqueNewMets)) = uniqueNewMets;
0471         [~, J]=ismember(newMets,allMets);
0472     end
0473 
0474     %Update the stoichiometric matrix for the model to add
0475     newS=sparse(numel(model.mets),numel(models{i}.rxns));
0476     newS(J,:)=models{i}.S;
0477     model.S=[model.S newS];
0478 
0479     
0480     %Now add new genes
0481     if isfield(models{i},'genes')
0482         if ~isfield(model,'genes')
0483             %If there was no gene info before
0484             model.genes=models{i}.genes;
0485             model.rxnGeneMat=[sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
0486             emptyGene=cell(numel(model.rxns),1);
0487             emptyGene(:)={''};
0488             model.grRules=[emptyGene;models{i}.grRules];
0489             model.geneFrom=cell(numel(models{i}.genes),1);
0490             model.geneFrom(:)={models{i}.id};
0491             
0492             if isfield(models{i},'geneShortNames')
0493                 model.geneShortNames=models{i}.geneShortNames;
0494             end
0495 
0496             if isfield(models{i},'proteins')
0497                 model.proteins=models{i}.proteins;
0498             end
0499 
0500             if isfield(models{i},'geneMiriams')
0501                 model.geneMiriams=models{i}.geneMiriams;
0502             end
0503             
0504             if isfield(models{i},'geneComps')
0505                 model.geneComps=models{i}.geneComps;
0506             end
0507         else
0508             %If gene info should be merged
0509             a=ismember(models{i}.genes,model.genes);
0510             
0511             genesToAdd=find(~a);
0512             
0513             %Only add extra gene info on new genes. This might not be
0514             %correct and should be changed later...
0515             if ~isempty(genesToAdd)
0516                 model.genes=[model.genes;models{i}.genes(genesToAdd)];
0517                 emptyGene=cell(numel(genesToAdd),1);
0518                 emptyGene(:)={models{i}.id};
0519                 model.geneFrom=[model.geneFrom;emptyGene];
0520                 model.rxnGeneMat=[model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
0521                 
0522                 if isfield(models{i},'geneShortNames')
0523                     if isfield(model,'geneShortNames')
0524                         model.geneShortNames=[model.geneShortNames;models{i}.geneShortNames(genesToAdd)];
0525                     else
0526                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
0527                         emptyGeneSN(:)={''};
0528                         model.geneShortNames=[emptyGeneSN;models{i}.geneShortNames(genesToAdd)];
0529                     end
0530                 else
0531                     if isfield(model,'geneShortNames')
0532                         emptyGeneSN=cell(numel(genesToAdd),1);
0533                         emptyGeneSN(:)={''};
0534                         model.geneShortNames=[model.geneShortNames;emptyGeneSN];
0535                     end
0536                 end
0537 
0538                 if isfield(models{i},'proteins')
0539                     if isfield(model,'proteins')
0540                         model.proteins=[model.proteins;models{i}.proteins(genesToAdd)];
0541                     else
0542                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
0543                         emptyGeneSN(:)={''};
0544                         model.proteins=[emptyGeneSN;models{i}.proteins(genesToAdd)];
0545                     end
0546                 else
0547                     if isfield(model,'proteins')
0548                         emptyGeneSN=cell(numel(genesToAdd),1);
0549                         emptyGeneSN(:)={''};
0550                         model.proteins=[model.proteins;emptyGeneSN];
0551                     end
0552                 end
0553 
0554                 if isfield(models{i},'geneMiriams')
0555                     if isfield(model,'geneMiriams')
0556                         model.geneMiriams=[model.geneMiriams;models{i}.geneMiriams(genesToAdd)];
0557                     else
0558                         emptyGeneMir=cell(numel(model.genes)-numel(genesToAdd),1);
0559                         model.geneMiriams=[emptyGeneMir;models{i}.geneMiriams(genesToAdd)];
0560                     end
0561                 else
0562                     if isfield(model,'geneMiriams')
0563                         emptyGeneMir=cell(numel(genesToAdd),1);
0564                         model.geneMiriams=[model.geneMiriams;emptyGeneMir];
0565                     end
0566                 end
0567                 
0568                 if isfield(models{i},'geneComps')
0569                     if isfield(model,'geneComps')
0570                         model.geneComps=[model.geneComps;models{i}.geneComps(genesToAdd)];
0571                     else
0572                         emptyGeneMir=ones(numel(model.genes)-numel(genesToAdd),1);
0573                         model.geneComps=[emptyGeneMir;models{i}.geneComps(genesToAdd)];
0574                         EM='Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment';
0575                         dispEM(EM,false);
0576                     end
0577                 else
0578                     if isfield(model,'geneComps')
0579                         emptyGeneMir=ones(numel(genesToAdd),1);
0580                         model.geneComps=[model.geneComps;emptyGeneMir];
0581                         EM='Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment';
0582                         dispEM(EM,false);
0583                     end
0584                 end
0585             end
0586             
0587             %Remap the genes from the new model. The same thing as with
0588             %mets; this is a wasteful way to do it but I don't care right
0589             %now
0590             [a, b]=ismember(models{i}.genes,model.genes);
0591             
0592             %Just a check
0593             if ~all(a)
0594                 EM='There was an unexpected error in matching genes';
0595                 dispEM(EM);
0596             end
0597             model.grRules=[model.grRules;models{i}.grRules];
0598         end
0599     else
0600         %Add empty gene associations
0601         if isfield(model,'genes')
0602             emptyGene=cell(numel(models{i}.rxns),1);
0603             emptyGene(:)={''};
0604             model.grRules=[model.grRules;emptyGene];
0605         end
0606     end
0607 end
0608 %Fix grRules and reconstruct rxnGeneMat
0609 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
0610 model.grRules = grRules;
0611 model.rxnGeneMat = rxnGeneMat;
0612 end