Home > RAVEN > mergeModels.m

mergeModels

PURPOSE ^

mergeModels

SYNOPSIS ^

function model=mergeModels(models,supressWarnings)

DESCRIPTION ^

 mergeModels
   Merges models into one model structure.

   models          a cell array with model structures
   supressWarnings true if warnings should be supressed (opt, 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)

   Rasmus Agren, 2013-02-07

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function model=mergeModels(models,supressWarnings)
0002 % mergeModels
0003 %   Merges models into one model structure.
0004 %
0005 %   models          a cell array with model structures
0006 %   supressWarnings true if warnings should be supressed (opt, default
0007 %                   false)
0008 %
0009 %   model     a model structure with the merged model. Follows the structure
0010 %             of normal models but also has 'rxnFrom/metFrom/geneFrom' fields
0011 %             to indicate from which model each reaction/metabolite/gene was
0012 %             taken
0013 %
0014 %   Usage: model=mergeModels(models)
0015 %
0016 %   Rasmus Agren, 2013-02-07
0017 %
0018 
0019 %Just return the model
0020 if numel(models)<=1
0021     model=models{1};
0022     return;
0023 end
0024 
0025 if nargin<2
0026     supressWarnings=false;
0027 end
0028 
0029 %Add new functionality in the order specified in models
0030 model=models{1};
0031 model.id='MERGED';
0032 model.description='';
0033 
0034 model.rxnFrom=cell(numel(models{1}.rxns),1);
0035 model.rxnFrom(:)={models{1}.id};
0036 model.metFrom=cell(numel(models{1}.mets),1);
0037 model.metFrom(:)={models{1}.id};
0038 if isfield(models{1},'genes')
0039     model.geneFrom=cell(numel(models{1}.genes),1);
0040     model.geneFrom(:)={models{1}.id};
0041 end
0042 
0043 if isfield(model,'subSystems')
0044     hasDeletedSubSystem=false;
0045 else
0046     if supressWarnings==false
0047         fprintf('WARNING: Cannot add subsystems since the existing model has no subsystems info. All reactions must have a subsystem for this to be included.\n');  
0048     end
0049     hasDeletedSubSystem=true;
0050 end
0051 
0052 %NOTE: I don't do this for rxnMiriam because it's not used at the moment
0053 if isfield(model,'rxnMiriams')
0054     model=rmfield(model,'rxnMiriams');
0055 end
0056 if isfield(model,'equations')
0057     model=rmfield(model,'equations');
0058 end
0059 
0060 for i=2:numel(models)
0061     %Add the model id to the rxn id id it already exists in the model (id
0062     %have to be unique)
0063     %This is because it makes a '[]' string if no new reactions
0064     if ~isempty(models{i}.rxns)
0065         I=ismember(models{i}.rxns,model.rxns);
0066         models{i}.rxns(I)=strcat(models{i}.rxns(I),['_' models{i}.id]);
0067     end 
0068     
0069     %Make sure that there are no conflicting reaction ids
0070     [crap crap conflicting]=intersect(model.rxns,models{i}.rxns);
0071     
0072     if ~isempty(conflicting)
0073        if supressWarnings==false
0074             fprintf(['WARNING: The following reaction IDs in ' models{i}.id ' are already present in the model and were renamed\n']);
0075        end
0076        for j=1:numel(conflicting)
0077            if supressWarnings==false
0078                 fprintf(['Old: ' models{i}.rxns{conflicting(j)} ' New: ' models{i}.rxns{conflicting(j)} '_' models{i}.id '\n']);
0079            end
0080            models{i}.rxns{conflicting(j)}=[models{i}.rxns{conflicting(j)} '_' models{i}.id];
0081        end
0082        if supressWarnings==false
0083             fprintf('\n');
0084        end
0085     end
0086     
0087     %Add all static stuff
0088     rxnFrom=cell(numel(models{i}.rxns),1);
0089     rxnFrom(:)={models{i}.id};
0090     model.rxnFrom=[model.rxnFrom;rxnFrom];
0091     model.rxns=[model.rxns;models{i}.rxns];
0092     model.rxnNames=[model.rxnNames;models{i}.rxnNames];
0093     model.lb=[model.lb;models{i}.lb];
0094     model.ub=[model.ub;models{i}.ub];
0095     model.c=[model.c;models{i}.c];
0096     model.rev=[model.rev;models{i}.rev];
0097     
0098     if hasDeletedSubSystem==false
0099         if isfield(models{i},'subSystems')
0100             model.subSystems=[model.subSystems;models{i}.subSystems];
0101         else
0102            if supressWarnings==false
0103                 fprintf('WARNING: Cannot add subsystems since the existing model has no subsystems info. All reactions must have a subsystem for this to be included. Deleting subSystems field.\n');  
0104            end
0105            hasDeletedSubSystem=true;
0106            model=rmfield(model,'subSystems');
0107         end
0108     end
0109     
0110     if isfield(models{i},'eccodes')
0111        if isfield(model,'eccodes')
0112            model.eccodes=[model.eccodes;models{i}.eccodes];
0113        else
0114            emptyEC=cell(numel(model.rxns)-numel(models{i}.rxns),1);
0115            emptyEC(:)={''};
0116            model.eccodes=[emptyEC;models{i}.eccodes];
0117        end
0118     else
0119        if isfield(model,'eccodes')
0120            emptyEC=cell(numel(models{i}.rxns),1);
0121            emptyEC(:)={''};
0122            model.eccodes=[model.eccodes;emptyEC];
0123        end 
0124     end
0125     
0126     if isfield(models{i},'rxnScores')
0127        if isfield(model,'rxnScores')
0128            model.rxnScores=[model.rxnScores;models{i}.rxnScores];
0129        else
0130            emptyRS=zeros(numel(model.rxns)-numel(models{i}.rxns),1);
0131            model.rxnScores=[emptyRS;models{i}.rxnScores];
0132        end
0133     else
0134        if isfield(model,'rxnScores')
0135            emptyRS=zeros(numel(models{i}.rxns),1);
0136            model.rxnScores=[model.rxnScores;emptyRS];
0137        end 
0138     end
0139     
0140     %Get the new metabolites from matching the models.
0141     %Metabolites are said to be the same if they share name and
0142     %compartment id. This means that metabolite IDs are not taken into
0143     %account.
0144     oldMetComps=model.comps(model.metComps);
0145     oldMets=strcat(model.metNames,'[',oldMetComps,']');
0146     %This is because it makes a '[]' string if no new metabolites
0147     if ~isempty(models{i}.metNames)
0148         newMetComps=models{i}.comps(models{i}.metComps);
0149         newMets=strcat(models{i}.metNames,'[',newMetComps,']');
0150     else
0151         newMets={};
0152         newMetComps={};
0153     end
0154     [tf,loc]=ismember(newMets,oldMets);
0155     metsToAdd=find(~tf);
0156     
0157     %First add the new metabolites
0158     %Make sure that there are no conflicting metabolite ids
0159     [conflicting,crap]=ismember(models{i}.mets(metsToAdd),model.mets);
0160     
0161     conflicting=find(conflicting);
0162     
0163     if ~isempty(conflicting)
0164        if supressWarnings==false 
0165             fprintf(['WARNING: The following metabolite IDs in ' models{i}.id ' are already present in the model and were renamed\n']);
0166        end
0167        for j=1:numel(conflicting)
0168            if supressWarnings==false
0169                 fprintf(['Old: ' models{i}.mets{metsToAdd(conflicting(j))} ' New: ' models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id '\n']);
0170            end
0171            models{i}.mets{metsToAdd(conflicting(j))}=[models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
0172        end
0173        if supressWarnings==false
0174             fprintf('\n');
0175        end
0176     end
0177     
0178     %Add static info on the metabolites
0179     metFrom=cell(numel(metsToAdd),1);
0180     metFrom(:)={models{i}.id};
0181     model.metFrom=[model.metFrom;metFrom];
0182     model.mets=[model.mets;models{i}.mets(metsToAdd)];
0183     model.metNames=[model.metNames;models{i}.metNames(metsToAdd)];
0184     model.b=[model.b;zeros(numel(metsToAdd),size(model.b,2))];
0185     
0186     if isfield(model,'unconstrained')
0187        if isfield(models{i},'unconstrained')
0188             model.unconstrained=[model.unconstrained;models{i}.unconstrained(metsToAdd)];
0189        else
0190            model.unconstrained=[model.unconstrained;zeros(numel(metsToAdd),1)];
0191        end
0192     else
0193        if isfield(models{i},'unconstrained')
0194           model.unconstrained=[zeros(numel(model.mets),1);models{i}.unconstrained(metsToAdd)];
0195        end
0196     end
0197     
0198     %Only add extra info on new metabolites since it's a little tricky to
0199     %chose what to keep otherwise. Should change in the future
0200     if ~isempty(metsToAdd)
0201         if isfield(models{i},'inchis')
0202            if isfield(model,'inchis')
0203                model.inchis=[model.inchis;models{i}.inchis(metsToAdd)];
0204            else
0205                emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
0206                emptyInchi(:)={''};
0207                model.inchis=[emptyInchi;models{i}.inchis(metsToAdd)];
0208            end
0209         else
0210            if isfield(model,'inchis')
0211                emptyInchi=cell(numel(metsToAdd),1);
0212                emptyInchi(:)={''};
0213                model.inchis=[model.inchis;emptyInchi];
0214            end 
0215         end
0216         
0217         if isfield(models{i},'metFormulas')
0218            if isfield(model,'metFormulas')
0219                model.metFormulas=[model.metFormulas;models{i}.metFormulas(metsToAdd)];
0220            else
0221                emptyMetFormulas=cell(numel(model.mets)-numel(metsToAdd),1);
0222                emptyMetFormulas(:)={''};
0223                model.metFormulas=[emptyMetFormulas;models{i}.metFormulas(metsToAdd)];
0224            end
0225         else
0226            if isfield(model,'metFormulas')
0227                emptyMetFormulas=cell(numel(metsToAdd),1);
0228                emptyMetFormulas(:)={''};
0229                model.metFormulas=[model.metFormulas;emptyMetFormulas];
0230            end 
0231         end
0232         
0233         if isfield(models{i},'metMiriams')
0234            if isfield(model,'metMiriams')
0235                model.metMiriams=[model.metMiriams;models{i}.metMiriams(metsToAdd)];
0236            else
0237                emptyMetMiriam=cell(numel(model.mets)-numel(metsToAdd),1);
0238                model.metMiriams=[emptyMetMiriam;models{i}.metMiriams(metsToAdd)];
0239            end
0240         else
0241            if isfield(model,'metMiriams')
0242                emptyMetMiriam=cell(numel(metsToAdd),1);
0243                model.metMiriams=[model.metMiriams;emptyMetMiriam];
0244            end 
0245         end
0246     end
0247     
0248     %Add if there are any new compartments and add those. This can change
0249     %the order of compartments and the corresponding indexes in
0250     %model.metComps.
0251     
0252     %Find overlapping and new compartments
0253     [overlap oldIDs]=ismember(models{i}.comps,model.comps);
0254     overlap=find(overlap);
0255     
0256     %Add the new compartments if any
0257     if numel(overlap)~=numel(models{i}.compNames)
0258         compIndexes=oldIDs==0;
0259         
0260         %Make sure that there are no conflicting compartment ids
0261         [conflicting, conflictID]=ismember(models{i}.compNames(compIndexes),model.compNames);
0262 
0263         conflicting=find(conflicting);
0264 
0265         if ~isempty(conflicting)
0266             message=['ERROR: The following compartment IDs in ' models{i}.id ' are already present in the model but with another name. They have to be renamed\n'];
0267             for j=1:numel(conflicting)
0268                 message=[message models{i}.comps{conflictID(conflicting(j))} '\n'];
0269             end
0270             throw(MException('',message));
0271         end
0272         
0273         %It's ok to add duplicate name, but not duplicate IDs
0274         model.compNames=[model.compNames; models{i}.compNames(compIndexes)];
0275         model.comps=[model.comps; models{i}.comps(compIndexes)];
0276         if isfield(model,'compOutside')
0277             model.compOutside=[model.compOutside; models{i}.compOutside(compIndexes)];
0278         end
0279     end
0280     
0281     %Only add new comp info on the un-matched metabolites since the old ones will
0282     %be mapped to the existing list anyways
0283     [I J]=ismember(newMetComps(metsToAdd),model.comps);
0284     %Just a check
0285     if ~all(I)
0286         throw(MException('','There was an unexplained error in matching compartments'));
0287     end
0288     model.metComps=[model.metComps;J];
0289     
0290     %Create the new stoichiometric matrix
0291     model.S=[model.S;sparse(numel(metsToAdd),size(model.S,2))];
0292     
0293     %Rematch metabolite names. Not the most clever way to do it maybe
0294     allMets=strcat(model.metNames,'[',model.comps(model.metComps),']');
0295     [I J]=ismember(newMets,allMets);
0296     
0297     %Update the stoichiometric matrix for the model to add
0298     newS=sparse(numel(model.mets),numel(models{i}.rxns));
0299     newS(J,:)=models{i}.S;
0300     
0301     model.S=[model.S newS];
0302     
0303     %Now add new genes
0304     if isfield(models{i},'genes')
0305         if ~isfield(model,'genes')
0306             %If there was no gene info before
0307             model.genes=models{i}.genes;
0308             model.rxnGeneMat=[sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
0309             emptyGene=cell(numel(model.rxns),1);
0310             emptyGene(:)={''};
0311             model.grRules=[emptyGene;models{i}.grRules];
0312             model.geneFrom=cell(numel(models{i}.genes),1);
0313             model.geneFrom(:)={models{i}.id};
0314             
0315             if isfield(models{i},'geneShortNames')
0316                model.geneShortNames=models{i}.geneShortNames; 
0317             end
0318             
0319             if isfield(models{i},'geneMiriams')
0320                model.geneMiriams=models{i}.geneMiriams; 
0321             end
0322         else
0323             %If gene info should be merged
0324             [a crap]=ismember(models{i}.genes,model.genes);
0325             
0326             genesToAdd=find(~a);
0327             
0328             %Only add extra gene info on new genes. This might not be
0329             %correct and should be changed later...
0330             if ~isempty(genesToAdd)
0331                 model.genes=[model.genes;models{i}.genes(genesToAdd)];
0332                 emptyGene=cell(numel(genesToAdd),1);
0333                 emptyGene(:)={models{i}.id};
0334                 model.geneFrom=[model.geneFrom;emptyGene];
0335                 model.rxnGeneMat=[model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
0336                 
0337                 if isfield(models{i},'geneShortNames')
0338                     if isfield(model,'geneShortNames')
0339                         model.geneShortNames=[model.geneShortNames;models{i}.geneShortNames(genesToAdd)];
0340                     else
0341                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
0342                         emptyGeneSN(:)={''};
0343                         model.geneShortNames=[emptyGeneSN;models{i}.geneShortNames(genesToAdd)];
0344                     end
0345                 else
0346                     if isfield(model,'geneShortNames')
0347                         emptyGeneSN=cell(numel(genesToAdd),1);
0348                         emptyGeneSN(:)={''};
0349                         model.geneShortNames=[model.geneShortNames;emptyGeneSN];   
0350                     end
0351                 end
0352 
0353                 if isfield(models{i},'geneMiriams')
0354                     if isfield(model,'geneMiriams')
0355                         model.geneMiriams=[model.geneMiriams;models{i}.geneMiriams(genesToAdd)];
0356                     else
0357                         emptyGeneMir=cell(numel(model.genes)-numel(genesToAdd),1);
0358                         model.geneMiriams=[emptyGeneMir;models{i}.geneMiriams(genesToAdd)];
0359                     end
0360                 else
0361                     if isfield(model,'geneMiriams')
0362                         emptyGeneMir=cell(numel(genesToAdd),1);
0363                         model.geneMiriams=[model.geneMiriams;emptyGeneMir];   
0364                     end
0365                 end
0366             end
0367             
0368             %Remap the genes from the new model. The same thing as with
0369             %mets; this is a wasteful way to do it but I don't care right
0370             %now
0371             [a b]=ismember(models{i}.genes,model.genes);
0372             
0373             %Just a check
0374             if ~all(a)
0375                 throw(MException('','There was an unexplained error in matching genes'));
0376             end
0377             
0378             %Create the new rxnGene matrix
0379             rxnGeneMat=sparse(numel(models{i}.rxns),numel(model.genes));
0380             rxnGeneMat(:,b)=models{i}.rxnGeneMat;
0381             model.rxnGeneMat=[model.rxnGeneMat; rxnGeneMat];
0382             model.grRules=[model.grRules;models{i}.grRules];
0383         end
0384     else
0385         %Add empty gene associations
0386         if isfield(model,'genes')
0387             model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(models{i}.rxns),numel(model.genes))];
0388             emptyGene=cell(numel(models{i}.rxns),1);
0389             emptyGene(:)={''};
0390             model.grRules=[model.grRules;emptyGene];
0391         end
0392     end
0393 end
0394 end

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005