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-08-01

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-08-01
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         dispEM('Cannot add subsystems since the existing model has no subsystems info. All reactions must have a subsystem for this to be included',false);  
0048     end
0049     hasDeletedSubSystem=true;
0050 end
0051 
0052 if isfield(model,'equations')
0053     model=rmfield(model,'equations');
0054 end
0055 
0056 for i=2:numel(models)
0057     %Add the model id to the rxn id id it already exists in the model (id
0058     %have to be unique)
0059     %This is because it makes a '[]' string if no new reactions
0060     if ~isempty(models{i}.rxns)
0061         I=ismember(models{i}.rxns,model.rxns);
0062         models{i}.rxns(I)=strcat(models{i}.rxns(I),['_' models{i}.id]);
0063     end 
0064     
0065     %Make sure that there are no conflicting reaction ids
0066     [crap crap conflicting]=intersect(model.rxns,models{i}.rxns);
0067     
0068     if ~isempty(conflicting)
0069        printString=cell(numel(conflicting),1);
0070        for j=1:numel(conflicting)
0071            printString{j}=['Old: ' models{i}.rxns{conflicting(j)} ' New: ' models{i}.rxns{conflicting(j)} '_' models{i}.id];
0072            models{i}.rxns{conflicting(j)}=[models{i}.rxns{conflicting(j)} '_' models{i}.id];
0073        end
0074        if supressWarnings==false
0075             dispEM(['The following reaction IDs in ' models{i}.id ' are already present in the model and were renamed:'],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 hasDeletedSubSystem==false
0092         if isfield(models{i},'subSystems')
0093             model.subSystems=[model.subSystems;models{i}.subSystems];
0094         else
0095            if supressWarnings==false
0096                 dispEM('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',false);  
0097            end
0098            hasDeletedSubSystem=true;
0099            model=rmfield(model,'subSystems');
0100         end
0101     end
0102     
0103     if isfield(models{i},'eccodes')
0104        if isfield(model,'eccodes')
0105            model.eccodes=[model.eccodes;models{i}.eccodes];
0106        else
0107            emptyEC=cell(numel(model.rxns)-numel(models{i}.rxns),1);
0108            emptyEC(:)={''};
0109            model.eccodes=[emptyEC;models{i}.eccodes];
0110        end
0111     else
0112        if isfield(model,'eccodes')
0113            emptyEC=cell(numel(models{i}.rxns),1);
0114            emptyEC(:)={''};
0115            model.eccodes=[model.eccodes;emptyEC];
0116        end 
0117     end
0118     
0119     if isfield(models{i},'rxnMiriams')
0120        if isfield(model,'rxnMiriams')
0121            model.rxnMiriams=[model.rxnMiriams;models{i}.rxnMiriams];
0122        else
0123            model.rxnMiriams=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnMiriams];
0124        end
0125     else
0126        if isfield(model,'rxnMiriams')
0127            model.rxnMiriams=[model.rxnMiriams;cell(numel(models{i}.rxns),1)];
0128        end 
0129     end
0130     
0131     if isfield(models{i},'rxnComps')
0132        if isfield(model,'rxnComps')
0133            model.rxnComps=[model.rxnComps;models{i}.rxnComps];
0134        else
0135            model.rxnComps=[ones(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnComps];
0136            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');
0137        end
0138     else
0139        if isfield(model,'rxnComps')
0140            model.rxnComps=[model.rxnComps;ones(numel(models{i}.rxns),1)];
0141            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');
0142        end 
0143     end
0144     
0145     if isfield(models{i},'rxnScores')
0146        if isfield(model,'rxnScores')
0147            model.rxnScores=[model.rxnScores;models{i}.rxnScores];
0148        else
0149            emptyRS=zeros(numel(model.rxns)-numel(models{i}.rxns),1);
0150            model.rxnScores=[emptyRS;models{i}.rxnScores];
0151        end
0152     else
0153        if isfield(model,'rxnScores')
0154            emptyRS=zeros(numel(models{i}.rxns),1);
0155            model.rxnScores=[model.rxnScores;emptyRS];
0156        end 
0157     end
0158     
0159     %Get the new metabolites from matching the models.
0160     %Metabolites are said to be the same if they share name and
0161     %compartment id. This means that metabolite IDs are not taken into
0162     %account.
0163     oldMetComps=model.comps(model.metComps);
0164     oldMets=strcat(model.metNames,'[',oldMetComps,']');
0165     %This is because it makes a '[]' string if no new metabolites
0166     if ~isempty(models{i}.metNames)
0167         newMetComps=models{i}.comps(models{i}.metComps);
0168         newMets=strcat(models{i}.metNames,'[',newMetComps,']');
0169     else
0170         newMets={};
0171         newMetComps={};
0172     end
0173     tf=ismember(newMets,oldMets);
0174     metsToAdd=find(~tf);
0175     
0176     %First add the new metabolites
0177     %Make sure that there are no conflicting metabolite ids
0178     [conflicting,crap]=ismember(models{i}.mets(metsToAdd),model.mets);
0179     
0180     conflicting=find(conflicting);
0181     
0182     if ~isempty(conflicting)
0183        printString=cell(numel(conflicting),1);
0184        for j=1:numel(conflicting)
0185            printString{j}=['Old: ' models{i}.mets{metsToAdd(conflicting(j))} ' New: ' models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
0186            models{i}.mets{metsToAdd(conflicting(j))}=[models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
0187        end
0188        if supressWarnings==false
0189            dispEM(['The following metabolite IDs in ' models{i}.id ' are already present in the model and were renamed:'],false,printString);
0190        end
0191     end
0192     
0193     %Add static info on the metabolites
0194     metFrom=cell(numel(metsToAdd),1);
0195     metFrom(:)={models{i}.id};
0196     model.metFrom=[model.metFrom;metFrom];
0197     model.mets=[model.mets;models{i}.mets(metsToAdd)];
0198     model.metNames=[model.metNames;models{i}.metNames(metsToAdd)];
0199     model.b=[model.b;zeros(numel(metsToAdd),size(model.b,2))];
0200     
0201     if isfield(model,'unconstrained')
0202        if isfield(models{i},'unconstrained')
0203             model.unconstrained=[model.unconstrained;models{i}.unconstrained(metsToAdd)];
0204        else
0205            model.unconstrained=[model.unconstrained;zeros(numel(metsToAdd),1)];
0206        end
0207     else
0208        if isfield(models{i},'unconstrained')
0209           model.unconstrained=[zeros(numel(model.mets),1);models{i}.unconstrained(metsToAdd)];
0210        end
0211     end
0212     
0213     %Only add extra info on new metabolites since it's a little tricky to
0214     %chose what to keep otherwise. Should change in the future
0215     if ~isempty(metsToAdd)
0216         if isfield(models{i},'inchis')
0217            if isfield(model,'inchis')
0218                model.inchis=[model.inchis;models{i}.inchis(metsToAdd)];
0219            else
0220                emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
0221                emptyInchi(:)={''};
0222                model.inchis=[emptyInchi;models{i}.inchis(metsToAdd)];
0223            end
0224         else
0225            if isfield(model,'inchis')
0226                emptyInchi=cell(numel(metsToAdd),1);
0227                emptyInchi(:)={''};
0228                model.inchis=[model.inchis;emptyInchi];
0229            end 
0230         end
0231         
0232         if isfield(models{i},'metFormulas')
0233            if isfield(model,'metFormulas')
0234                model.metFormulas=[model.metFormulas;models{i}.metFormulas(metsToAdd)];
0235            else
0236                emptyMetFormulas=cell(numel(model.mets)-numel(metsToAdd),1);
0237                emptyMetFormulas(:)={''};
0238                model.metFormulas=[emptyMetFormulas;models{i}.metFormulas(metsToAdd)];
0239            end
0240         else
0241            if isfield(model,'metFormulas')
0242                emptyMetFormulas=cell(numel(metsToAdd),1);
0243                emptyMetFormulas(:)={''};
0244                model.metFormulas=[model.metFormulas;emptyMetFormulas];
0245            end 
0246         end
0247         
0248         if isfield(models{i},'metMiriams')
0249            if isfield(model,'metMiriams')
0250                model.metMiriams=[model.metMiriams;models{i}.metMiriams(metsToAdd)];
0251            else
0252                emptyMetMiriam=cell(numel(model.mets)-numel(metsToAdd),1);
0253                model.metMiriams=[emptyMetMiriam;models{i}.metMiriams(metsToAdd)];
0254            end
0255         else
0256            if isfield(model,'metMiriams')
0257                emptyMetMiriam=cell(numel(metsToAdd),1);
0258                model.metMiriams=[model.metMiriams;emptyMetMiriam];
0259            end 
0260         end
0261     end
0262     
0263     %Add if there are any new compartments and add those. This can change
0264     %the order of compartments and the corresponding indexes in
0265     %model.metComps.
0266     
0267     %Find overlapping and new compartments
0268     [overlap oldIDs]=ismember(models{i}.comps,model.comps);
0269     overlap=find(overlap);
0270     
0271     %Add the new compartments if any
0272     if numel(overlap)~=numel(models{i}.compNames)
0273         compIndexes=oldIDs==0;
0274         
0275         %Make sure that there are no conflicting compartment ids
0276         [crap, conflicting]=ismember(models{i}.compNames(compIndexes),model.compNames);
0277         dispEM(['The following compartment IDs in ' models{i}.id ' are already present in the model but with another name. They have to be renamed'],true,model.comps(conflicting));
0278         
0279         %It's ok to add duplicate name, but not duplicate IDs
0280         model.compNames=[model.compNames; models{i}.compNames(compIndexes)];
0281         model.comps=[model.comps; models{i}.comps(compIndexes)];
0282         if isfield(model,'compOutside')
0283             if isfield(models{i},'compOutside')
0284                 model.compOutside=[model.compOutside; models{i}.compOutside(compIndexes)];
0285             else
0286                 %This is if not all models have the field
0287                 padding=cell(sum(compIndexes),1);
0288                 padding(:)={''};
0289                 model.compOutside=[model.compOutside;padding];
0290             end
0291         end
0292         if isfield(model,'compMiriams')
0293             if isfield(models{i},'compMiriams')
0294                 model.compMiriams=[model.compMiriams; models{i}.compMiriams(compIndexes)];
0295             else
0296                 %This is if not all models have the field
0297                 model.compMiriams=[model.compMiriams;cell(sum(compIndexes),1)];
0298             end
0299         end
0300     end
0301     
0302     %Only add new comp info on the un-matched metabolites since the old ones will
0303     %be mapped to the existing list anyways
0304     [I J]=ismember(newMetComps(metsToAdd),model.comps);
0305     %Just a check
0306     if ~all(I)
0307         dispEM('There was an unexplained error in matching compartments');
0308     end
0309     model.metComps=[model.metComps;J];
0310     
0311     %Create the new stoichiometric matrix
0312     model.S=[model.S;sparse(numel(metsToAdd),size(model.S,2))];
0313     
0314     %Rematch metabolite names. Not the most clever way to do it maybe
0315     allMets=strcat(model.metNames,'[',model.comps(model.metComps),']');
0316     [crap J]=ismember(newMets,allMets);
0317     
0318     %Update the stoichiometric matrix for the model to add
0319     newS=sparse(numel(model.mets),numel(models{i}.rxns));
0320     newS(J,:)=models{i}.S;
0321     
0322     model.S=[model.S newS];
0323     
0324     %Now add new genes
0325     if isfield(models{i},'genes')
0326         if ~isfield(model,'genes')
0327             %If there was no gene info before
0328             model.genes=models{i}.genes;
0329             model.rxnGeneMat=[sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
0330             emptyGene=cell(numel(model.rxns),1);
0331             emptyGene(:)={''};
0332             model.grRules=[emptyGene;models{i}.grRules];
0333             model.geneFrom=cell(numel(models{i}.genes),1);
0334             model.geneFrom(:)={models{i}.id};
0335             
0336             if isfield(models{i},'geneShortNames')
0337                model.geneShortNames=models{i}.geneShortNames; 
0338             end
0339             
0340             if isfield(models{i},'geneMiriams')
0341                model.geneMiriams=models{i}.geneMiriams; 
0342             end
0343             
0344             if isfield(models{i},'geneComps')
0345                model.geneComps=models{i}.geneComps; 
0346             end
0347         else
0348             %If gene info should be merged
0349             [a crap]=ismember(models{i}.genes,model.genes);
0350             
0351             genesToAdd=find(~a);
0352             
0353             %Only add extra gene info on new genes. This might not be
0354             %correct and should be changed later...
0355             if ~isempty(genesToAdd)
0356                 model.genes=[model.genes;models{i}.genes(genesToAdd)];
0357                 emptyGene=cell(numel(genesToAdd),1);
0358                 emptyGene(:)={models{i}.id};
0359                 model.geneFrom=[model.geneFrom;emptyGene];
0360                 model.rxnGeneMat=[model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
0361                 
0362                 if isfield(models{i},'geneShortNames')
0363                     if isfield(model,'geneShortNames')
0364                         model.geneShortNames=[model.geneShortNames;models{i}.geneShortNames(genesToAdd)];
0365                     else
0366                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
0367                         emptyGeneSN(:)={''};
0368                         model.geneShortNames=[emptyGeneSN;models{i}.geneShortNames(genesToAdd)];
0369                     end
0370                 else
0371                     if isfield(model,'geneShortNames')
0372                         emptyGeneSN=cell(numel(genesToAdd),1);
0373                         emptyGeneSN(:)={''};
0374                         model.geneShortNames=[model.geneShortNames;emptyGeneSN];   
0375                     end
0376                 end
0377 
0378                 if isfield(models{i},'geneMiriams')
0379                     if isfield(model,'geneMiriams')
0380                         model.geneMiriams=[model.geneMiriams;models{i}.geneMiriams(genesToAdd)];
0381                     else
0382                         emptyGeneMir=cell(numel(model.genes)-numel(genesToAdd),1);
0383                         model.geneMiriams=[emptyGeneMir;models{i}.geneMiriams(genesToAdd)];
0384                     end
0385                 else
0386                     if isfield(model,'geneMiriams')
0387                         emptyGeneMir=cell(numel(genesToAdd),1);
0388                         model.geneMiriams=[model.geneMiriams;emptyGeneMir];   
0389                     end
0390                 end
0391                 
0392                 if isfield(models{i},'geneComps')
0393                     if isfield(model,'geneComps')
0394                         model.geneComps=[model.geneComps;models{i}.geneComps(genesToAdd)];
0395                     else
0396                         emptyGeneMir=ones(numel(model.genes)-numel(genesToAdd),1);
0397                         model.geneComps=[emptyGeneMir;models{i}.geneComps(genesToAdd)];
0398                         dispEM('Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment',false);
0399                     end
0400                 else
0401                     if isfield(model,'geneComps')
0402                         emptyGeneMir=ones(numel(genesToAdd),1);
0403                         model.geneComps=[model.geneComps;emptyGeneMir];
0404                         dispEM('Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment',false);
0405                     end
0406                 end
0407             end
0408             
0409             %Remap the genes from the new model. The same thing as with
0410             %mets; this is a wasteful way to do it but I don't care right
0411             %now
0412             [a b]=ismember(models{i}.genes,model.genes);
0413             
0414             %Just a check
0415             if ~all(a)
0416                 dispEM('There was an unexplained error in matching genes');
0417             end
0418             
0419             %Create the new rxnGene matrix
0420             rxnGeneMat=sparse(numel(models{i}.rxns),numel(model.genes));
0421             rxnGeneMat(:,b)=models{i}.rxnGeneMat;
0422             model.rxnGeneMat=[model.rxnGeneMat; rxnGeneMat];
0423             model.grRules=[model.grRules;models{i}.grRules];
0424         end
0425     else
0426         %Add empty gene associations
0427         if isfield(model,'genes')
0428             model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(models{i}.rxns),numel(model.genes))];
0429             emptyGene=cell(numel(models{i}.rxns),1);
0430             emptyGene(:)={''};
0431             model.grRules=[model.grRules;emptyGene];
0432         end
0433     end
0434 end
0435 end

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005