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.

 Input:
   models          a cell array with model structures
   metParam        string metabolite name ('metNames') or ID ('mets') are
                   used for matching (optional, default 'metNames')
   supressWarnings logical whether warnings should be supressed (optional,
                   default false)
   copyToComps     logical whether mergeModels is run via copyToComps
                   (optional, default false)

 Output:
   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. If the model
                   already has 'rxnFrom/metFrom/geneFrom' fields, then
                   these fields are not modified.

 Usage: model=mergeModels(models)

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