From c322f335b9475ec7d0704e8fec02a90a3f4d4537 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 12 Nov 2025 23:27:26 +0100 Subject: [PATCH 1/9] fix: predictLocalization with CELLO output --- core/predictLocalization.m | 2 +- doc/core/predictLocalization.html | 2 +- doc/external/parseScores.html | 2 +- external/parseScores.m | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/core/predictLocalization.m b/core/predictLocalization.m index 2bfff833..a2b5a550 100755 --- a/core/predictLocalization.m +++ b/core/predictLocalization.m @@ -621,7 +621,7 @@ outModel.compNames(2)=GSS.compartments(1); end end -outModel.compNames=[outModel.compNames;GSS.compartments(2:end)']; +outModel.compNames=[outModel.compNames;GSS.compartments(2:end)]; %Ugly little loop for i=1:numel(GSS.compartments)-1 diff --git a/doc/core/predictLocalization.html b/doc/core/predictLocalization.html index 2622161e..317e1f4f 100644 --- a/doc/core/predictLocalization.html +++ b/doc/core/predictLocalization.html @@ -726,7 +726,7 @@

SOURCE CODE ^end 0623 end -0624 outModel.compNames=[outModel.compNames;GSS.compartments(2:end)']; +0624 outModel.compNames=[outModel.compNames;GSS.compartments(2:end)]; 0625 0626 %Ugly little loop 0627 for i=1:numel(GSS.compartments)-1 diff --git a/doc/external/parseScores.html b/doc/external/parseScores.html index 70766632..d82f859d 100644 --- a/doc/external/parseScores.html +++ b/doc/external/parseScores.html @@ -150,7 +150,7 @@

SOURCE CODE ^%Read the title line and fetch the list of compartments 0093 tline = fgetl(fid); 0094 GSS.compartments=regexp(tline,',','split'); -0095 GSS.compartments=GSS.compartments(4:end); +0095 GSS.compartments=transpose(GSS.compartments(4:end)); 0096 0097 %Now iterate through the following lines in the file. Each row 0098 %corresponds to one gene and it consists of the scores for diff --git a/external/parseScores.m b/external/parseScores.m index 662320c2..ae865b26 100755 --- a/external/parseScores.m +++ b/external/parseScores.m @@ -92,7 +92,7 @@ %Read the title line and fetch the list of compartments tline = fgetl(fid); GSS.compartments=regexp(tline,',','split'); - GSS.compartments=GSS.compartments(4:end); + GSS.compartments=transpose(GSS.compartments(4:end)); %Now iterate through the following lines in the file. Each row %corresponds to one gene and it consists of the scores for From dd6a816cb651f4383a245f116b11c5a38c317300 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 4 Dec 2025 15:46:13 +0100 Subject: [PATCH 2/9] fix: setParam bug --- core/setParam.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/setParam.m b/core/setParam.m index 734853b7..70705f3f 100755 --- a/core/setParam.m +++ b/core/setParam.m @@ -68,7 +68,7 @@ params(~Lia)=[]; indexes(~Lia)=[]; paramType(~Lia)=[]; - dispEM('Reactions not present in model, will be ignored:',false,rxnLise(~Lia)); + dispEM('Reactions not present in model, will be ignored:',false,rxnList(~Lia)); end %Change the parameters From 9e0c20aa2a67a07174c3c33e2ac7e863db020e84 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 4 Dec 2025 15:46:31 +0100 Subject: [PATCH 3/9] WIP: mergeModels copyToComps --- core/copyToComps.m | 14 ++++++++++---- core/mergeModels.m | 25 +++++++++++++++++-------- 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/core/copyToComps.m b/core/copyToComps.m index a285e650..84dad6bb 100755 --- a/core/copyToComps.m +++ b/core/copyToComps.m @@ -80,14 +80,20 @@ end modelToAdd.metComps=ones(numel(modelToAdd.mets),1); + if isfield(modelToAdd,'rxnFrom') + modelToAdd=rmfield(modelToAdd,'rxnFrom'); + end + if isfield(modelToAdd,'metFrom') + modelToAdd=rmfield(modelToAdd,'metFrom'); + end + if isfield(modelToAdd,'geneFrom') + modelToAdd=rmfield(modelToAdd,'geneFrom'); + end + %Merge the models model=mergeModels({model;modelToAdd},'metNames'); end -model=rmfield(model,'rxnFrom'); -model=rmfield(model,'metFrom'); -model=rmfield(model,'geneFrom'); - if deleteOriginal==true model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments end diff --git a/core/mergeModels.m b/core/mergeModels.m index 29faff96..586ca069 100755 --- a/core/mergeModels.m +++ b/core/mergeModels.m @@ -5,16 +5,20 @@ % their name and compartment (metaboliteName[comp]), while genes are % matched by their name. % +% Input: % 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) +% metParam string metabolite name ('metNames') or ID ('mets') are +% used for matching (default, metNames) +% supressWarnings logical whether 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 +% 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) @@ -34,11 +38,16 @@ supressWarnings=false; end +hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models);hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models); +hasGeneFrom = cellfun(@(s) isfield(s,'geneFrom'), models); +hasRxnsFrom = cellfun(@(s) isfield(s,'rxnsFrom'), models); + %Add new functionality in the order specified in models model=models{1}; model.id='MERGED'; model.name=''; +if ~any(hasRxnsFrom) model.rxnFrom=cell(numel(models{1}.rxns),1); model.rxnFrom(:)={models{1}.id}; model.metFrom=cell(numel(models{1}.mets),1); From e0eaa648781fe6239dd1007884271756f2680134 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 5 Dec 2025 22:22:16 +0100 Subject: [PATCH 4/9] fix: copyToComps keeps ...From fields --- core/copyToComps.m | 23 +- core/mergeModels.m | 80 +-- doc/core/copyToComps.html | 47 +- doc/core/mergeModels.html | 1227 +++++++++++++++++++------------------ doc/core/setParam.html | 2 +- 5 files changed, 708 insertions(+), 671 deletions(-) diff --git a/core/copyToComps.m b/core/copyToComps.m index 84dad6bb..37dd4ccb 100755 --- a/core/copyToComps.m +++ b/core/copyToComps.m @@ -79,21 +79,30 @@ modelToAdd.compMiriams=modelToAdd.compMiriams(J); end modelToAdd.metComps=ones(numel(modelToAdd.mets),1); - - if isfield(modelToAdd,'rxnFrom') - modelToAdd=rmfield(modelToAdd,'rxnFrom'); - end if isfield(modelToAdd,'metFrom') - modelToAdd=rmfield(modelToAdd,'metFrom'); + modelToAdd = rmfield(modelToAdd,'metFrom'); + end + if isfield(modelToAdd,'rxnFrom') + modelToAdd = rmfield(modelToAdd,'rxnFrom'); end if isfield(modelToAdd,'geneFrom') - modelToAdd=rmfield(modelToAdd,'geneFrom'); + modelToAdd = rmfield(modelToAdd,'geneFrom'); end - + %Merge the models model=mergeModels({model;modelToAdd},'metNames'); end +if all(strcmp(model.rxnFrom,originalID)) + model = rmfield(model,'rxnFrom'); +end +if all(strcmp(model.metFrom,originalID)) + model = rmfield(model,'metFrom'); +end +if all(strcmp(model.geneFrom,originalID)) + model = rmfield(model,'geneFrom'); +end + if deleteOriginal==true model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments end diff --git a/core/mergeModels.m b/core/mergeModels.m index 586ca069..47ad5007 100755 --- a/core/mergeModels.m +++ b/core/mergeModels.m @@ -38,25 +38,33 @@ supressWarnings=false; end -hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models);hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models); +hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models); hasGeneFrom = cellfun(@(s) isfield(s,'geneFrom'), models); -hasRxnsFrom = cellfun(@(s) isfield(s,'rxnsFrom'), models); +hasRxnFrom = cellfun(@(s) isfield(s,'rxnFrom'), models); + +for i = 1:numel(models) + if ~any(hasMetFrom) + models{i}.metFrom = repmat({models{i}.id},numel(models{i}.mets),1); + elseif ~hasMetFrom(i) + models{i}.metFrom = repmat({''},numel(models{i}.mets),1); + end + if ~any(hasRxnFrom) + models{i}.rxnFrom = repmat({models{i}.id},numel(models{i}.rxns),1); + elseif ~hasRxnFrom(i) + models{i}.rxnFrom = repmat({''},numel(models{i}.rxns),1); + end + if ~any(hasGeneFrom) && any(cellfun(@(s) isfield(s,'genes'), models)) + models{i}.geneFrom = repmat({models{i}.id},numel(models{i}.genes),1); + elseif ~hasGeneFrom(i) + models{i}.geneFrom = repmat({''},numel(models{i}.genes),1); + end +end %Add new functionality in the order specified in models model=models{1}; model.id='MERGED'; model.name=''; -if ~any(hasRxnsFrom) -model.rxnFrom=cell(numel(models{1}.rxns),1); -model.rxnFrom(:)={models{1}.id}; -model.metFrom=cell(numel(models{1}.mets),1); -model.metFrom(:)={models{1}.id}; -if isfield(models{1},'genes') - model.geneFrom=cell(numel(models{1}.genes),1); - model.geneFrom(:)={models{1}.id}; -end - if isfield(model,'equations') model=rmfield(model,'equations'); end @@ -87,15 +95,13 @@ end %Add all static stuff - rxnFrom=cell(numel(models{i}.rxns),1); - rxnFrom(:)={models{i}.id}; - model.rxnFrom=[model.rxnFrom;rxnFrom]; - model.rxns=[model.rxns;models{i}.rxns]; - model.rxnNames=[model.rxnNames;models{i}.rxnNames]; - model.lb=[model.lb;models{i}.lb]; - model.ub=[model.ub;models{i}.ub]; - model.c=[model.c;models{i}.c]; - model.rev=[model.rev;models{i}.rev]; + model.rxnFrom = [model.rxnFrom; models{i}.rxnFrom]; + model.rxns = [model.rxns; models{i}.rxns]; + model.rxnNames = [model.rxnNames; models{i}.rxnNames]; + model.lb = [model.lb; models{i}.lb]; + model.ub = [model.ub; models{i}.ub]; + model.c = [model.c; models{i}.c]; + model.rev = [model.rev; models{i}.rev]; if isfield(models{i},'subSystems') if isfield(model,'subSystems') @@ -296,12 +302,10 @@ end %Add static info on the metabolites - metFrom=cell(numel(metsToAdd),1); - metFrom(:)={models{i}.id}; - model.metFrom=[model.metFrom;metFrom]; - model.mets=[model.mets;models{i}.mets(metsToAdd)]; - model.metNames=[model.metNames;models{i}.metNames(metsToAdd)]; - model.b=[model.b;zeros(numel(metsToAdd),size(model.b,2))]; + model.metFrom = [model.metFrom; models{i}.metFrom(metsToAdd)]; + model.mets = [model.mets; models{i}.mets(metsToAdd)]; + model.metNames = [model.metNames; models{i}.metNames(metsToAdd)]; + model.b = [model.b; zeros(numel(metsToAdd),size(model.b,2))]; if isfield(model,'unconstrained') if isfield(models{i},'unconstrained') @@ -490,13 +494,11 @@ if isfield(models{i},'genes') if ~isfield(model,'genes') %If there was no gene info before - model.genes=models{i}.genes; - model.rxnGeneMat=[sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat]; - emptyGene=cell(numel(model.rxns),1); - emptyGene(:)={''}; - model.grRules=[emptyGene;models{i}.grRules]; - model.geneFrom=cell(numel(models{i}.genes),1); - model.geneFrom(:)={models{i}.id}; + model.genes = models{i}.genes; + model.rxnGeneMat = [sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat]; + emptyGene = repmat({''},numel(model.rxns),1); + model.grRules = [emptyGene;models{i}.grRules]; + model.geneFrom = models{i}.geneFrom; if isfield(models{i},'geneShortNames') model.geneShortNames=models{i}.geneShortNames; @@ -522,11 +524,9 @@ %Only add extra gene info on new genes. This might not be %correct and should be changed later... if ~isempty(genesToAdd) - model.genes=[model.genes;models{i}.genes(genesToAdd)]; - emptyGene=cell(numel(genesToAdd),1); - emptyGene(:)={models{i}.id}; - model.geneFrom=[model.geneFrom;emptyGene]; - model.rxnGeneMat=[model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))]; + model.genes = [model.genes; models{i}.genes(genesToAdd)]; + model.geneFrom = [model.geneFrom; models{i}.geneFrom(genesToAdd)]; + model.rxnGeneMat = [model.rxnGeneMat sparse(size(model.rxnGeneMat,1),numel(genesToAdd))]; if isfield(models{i},'geneShortNames') if isfield(model,'geneShortNames') @@ -596,7 +596,7 @@ %Remap the genes from the new model. The same thing as with %mets; this is a wasteful way to do it but I don't care right %now - [a, b]=ismember(models{i}.genes,model.genes); + a = ismember(models{i}.genes,model.genes); %Just a check if ~all(a) diff --git a/doc/core/copyToComps.html b/doc/core/copyToComps.html index 4e85433e..20c77753 100644 --- a/doc/core/copyToComps.html +++ b/doc/core/copyToComps.html @@ -147,22 +147,37 @@

SOURCE CODE ^end 0081 modelToAdd.metComps=ones(numel(modelToAdd.mets),1); -0082 -0083 %Merge the models -0084 model=mergeModels({model;modelToAdd},'metNames'); -0085 end -0086 -0087 model=rmfield(model,'rxnFrom'); -0088 model=rmfield(model,'metFrom'); -0089 model=rmfield(model,'geneFrom'); -0090 -0091 if deleteOriginal==true -0092 model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments -0093 end -0094 -0095 model.id=originalID; -0096 model.name=originalName; -0097 end +0082 if isfield(modelToAdd,'metFrom') +0083 modelToAdd = rmfield(modelToAdd,'metFrom'); +0084 end +0085 if isfield(modelToAdd,'rxnFrom') +0086 modelToAdd = rmfield(modelToAdd,'rxnFrom'); +0087 end +0088 if isfield(modelToAdd,'geneFrom') +0089 modelToAdd = rmfield(modelToAdd,'geneFrom'); +0090 end +0091 +0092 %Merge the models +0093 model=mergeModels({model;modelToAdd},'metNames'); +0094 end +0095 +0096 if all(strcmp(model.rxnFrom,originalID)) +0097 model = rmfield(model,'rxnFrom'); +0098 end +0099 if all(strcmp(model.metFrom,originalID)) +0100 model = rmfield(model,'metFrom'); +0101 end +0102 if all(strcmp(model.geneFrom,originalID)) +0103 model = rmfield(model,'geneFrom'); +0104 end +0105 +0106 if deleteOriginal==true +0107 model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments +0108 end +0109 +0110 model.id=originalID; +0111 model.name=originalName; +0112 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/core/mergeModels.html b/doc/core/mergeModels.html index b6d69f7c..530ecf62 100644 --- a/doc/core/mergeModels.html +++ b/doc/core/mergeModels.html @@ -33,16 +33,20 @@

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

SOURCE CODE ^dispEM('Reactions not present in model, will be ignored:',false,rxnLise(~Lia)); +0071 dispEM('Reactions not present in model, will be ignored:',false,rxnList(~Lia)); 0072 end 0073 0074 %Change the parameters From d79578a20a3512d4a6ff89989d3abe040def6411 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 22 Dec 2025 17:10:20 +0100 Subject: [PATCH 5/9] chore: fitTasks correct flags to fillGaps --- core/fitTasks.m | 2 +- core/parseTaskList.m | 9 +- doc/core/fitTasks.html | 2 +- doc/core/parseTaskList.html | 570 ++++++++++++++++++------------------ 4 files changed, 293 insertions(+), 290 deletions(-) diff --git a/core/fitTasks.m b/core/fitTasks.m index 9963f9e7..992c3762 100755 --- a/core/fitTasks.m +++ b/core/fitTasks.m @@ -266,7 +266,7 @@ %Only do gap-filling if it cannot be solved failed=false; try - [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params); + [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores); if exitFlag==-2 EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n']; dispEM(EM,false); diff --git a/core/parseTaskList.m b/core/parseTaskList.m index 614710c5..e6cafca0 100755 --- a/core/parseTaskList.m +++ b/core/parseTaskList.m @@ -2,10 +2,11 @@ % parseTaskList % Parses a task list file. % -% inputFile a task list in Excel format. The file must contain a -% sheet named TASKS, which in turn may contain the -% following column headers (note, all rows starting with -% a non-empty cell are removed. The first row after that +% inputFile a task list in either Excel (*.xlsx, with a sheet named +% TASKS with all relevant content) or tab-delimited +% (*.txt) format. The file may contain the following +% column headers (note, all rows starting with a +% non-empty cell are removed. The first row after that % is considered the headers): % ID % the only required header. Each task must have a diff --git a/doc/core/fitTasks.html b/doc/core/fitTasks.html index fa6bdbcb..fa51abdd 100644 --- a/doc/core/fitTasks.html +++ b/doc/core/fitTasks.html @@ -345,7 +345,7 @@

SOURCE CODE ^%Only do gap-filling if it cannot be solved 0267 failed=false; 0268 try -0269 [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores,params); +0269 [~, ~, newRxns, newModel, exitFlag]=fillGaps(tModel,refModel,false,true,supressWarnings,rxnScores); 0270 if exitFlag==-2 0271 EM=['"[' taskStructure(i).id '] ' taskStructure(i).description '" was aborted before reaching optimality. Consider increasing params.maxTime\n']; 0272 dispEM(EM,false); diff --git a/doc/core/parseTaskList.html b/doc/core/parseTaskList.html index 091bf07f..51aa28e2 100644 --- a/doc/core/parseTaskList.html +++ b/doc/core/parseTaskList.html @@ -30,10 +30,11 @@

DESCRIPTION ^
 parseTaskList
    Parses a task list file.
 
-   inputFile       a task list in Excel format. The file must contain a
-                   sheet named TASKS, which in turn may contain the
-                   following column headers (note, all rows starting with
-                   a non-empty cell are removed. The first row after that
+   inputFile       a task list in either Excel (*.xlsx, with a sheet named
+                   TASKS with all relevant content) or tab-delimited
+                   (*.txt) format. The file may contain the following
+                   column headers (note, all rows starting with a
+                   non-empty cell are removed. The first row after that
                    is considered the headers):
                    ID
                        the only required header. Each task must have a
@@ -154,286 +155,287 @@ 
SOURCE CODE ^% parseTaskList
 0003 %   Parses a task list file.
 0004 %
-0005 %   inputFile       a task list in Excel format. The file must contain a
-0006 %                   sheet named TASKS, which in turn may contain the
-0007 %                   following column headers (note, all rows starting with
-0008 %                   a non-empty cell are removed. The first row after that
-0009 %                   is considered the headers):
-0010 %                   ID
-0011 %                       the only required header. Each task must have a
-0012 %                       unique id (string or numeric). Tasks can span multiple
-0013 %                       rows, only the first row in each task should have
-0014 %                       an id
-0015 %                   DESCRIPTION
-0016 %                       description of the task
-0017 %                   IN
-0018 %                       allowed input(s) for the task. Metabolite names
-0019 %                       should be on the form
-0020 %                       "model.metName[model.comps]". Several inputs
-0021 %                       can be delimited by ";". If so, then the same
-0022 %                       bounds are used for all inputs. If that is not
-0023 %                       wanted, then use several rows for the task
-0024 %                   IN LB
-0025 %                       lower bound for the uptake of the metabolites in
-0026 %                       the row (optional, default 0 which corresponds to a
-0027 %                       minimal uptake of 0 units)
-0028 %                   IN UB
-0029 %                       upper bound for the uptake of the metabolites in
-0030 %                       the row (optional, default 1000 which corresponds to a
-0031 %                       maximal uptake of 1000 units)
-0032 %                   OUT
-0033 %                       allowed output(s) for the task (see IN)
-0034 %                   OUT LB
-0035 %                       lower bound for the production of the metabolites in
-0036 %                       the row (optional, default 0 which corresponds to a
-0037 %                       minimal production of 0 units)
-0038 %                   OUT UB
-0039 %                       upper bound for the production of the metabolites in
-0040 %                       the row (optional, default 1000 which corresponds to a
-0041 %                       maximal production of 1000 units)
-0042 %                   EQU
-0043 %                       equation to add. The equation should be on the form
-0044 %                       "0.4 A + 2 B <=> (or =>) C" and the metabolites
-0045 %                       should be on the form
-0046 %                       "model.metName[model.comps]" (optional)
-0047 %                   EQU LB
-0048 %                       lower bound for the equation (optional, default -1000
-0049 %                       for reversible and 0 for irreversible)
-0050 %                   EQU UB
-0051 %                       upper bound for the equation (optional, default 1000)
-0052 %                   CHANGED RXN
-0053 %                       reaction ID for which to change the bounds for.
-0054 %                       Several IDs can be delimited by ";". If so,
-0055 %                       then the same bounds are used for all reactions. If
-0056 %                       that is not wanted, then use several rows for the task
-0057 %                   CHANGED LB
-0058 %                       lower bound for the reaction
-0059 %                   CHANGED UB
-0060 %                       upper bound for the reaction
-0061 %                   SHOULD FAIL
-0062 %                       true if the correct behavior of the model is to
-0063 %                       not have a feasible solution given the constraints
-0064 %                       (optional, default false)
-0065 %                   PRINT FLUX
-0066 %                       true if the function should print the corresponding
-0067 %                       flux distribution for a task. Can be useful for
-0068 %                       testing (optional, default false)
-0069 %
-0070 %   taskStruct      array of structures with the following fields
-0071 %       id          the id of the task
-0072 %       description the description of the task
-0073 %       shouldFail  true if the task should fail
-0074 %       printFluxes true if the fluxes should be printed
-0075 %       comments    string with comments
-0076 %       inputs      cell array with input metabolites (in the form metName[comps])
-0077 %       LBin        array with lower bounds on inputs (default, 0)
-0078 %       UBin        array with upper bounds on inputs (default, 1000)
-0079 %       outputs     cell array with output metabolites (in the form metName[comps])
-0080 %       LBout       array with lower bounds on outputs (default, 0)
-0081 %       UBout       array with upper bounds on outputs (default, 1000)
-0082 %       equations   cell array with equations (with mets in the form metName[comps])
-0083 %       LBequ       array with lower bounds on equations (default, -1000 for
-0084 %                   reversible and 0 for irreversible)
-0085 %       UBequ       array with upper bounds on equations (default, 1000)
-0086 %       changed     cell array with reactions to change bounds for
-0087 %       LBrxn       array with lower bounds on changed reactions
-0088 %       UBrxn       array with upper bounds on changed reactions
-0089 %
-0090 %   This function is used for defining a set of tasks for a model to
-0091 %   perform. The tasks are defined by defining constraints on the model,
-0092 %   and if the problem is feasible, then the task is considered successful.
-0093 %   In general, each row can contain one constraint on uptakes, one
-0094 %   constraint on outputs, one new equation, and one change of reaction
-0095 %   bounds. If more bounds are needed to define the task, then several rows
-0096 %   can be used for each task. To perform the task use checkTasks or
-0097 %   fitTasks.
-0098 %
-0099 %   NOTE: The general metabolites "ALLMETS" and "ALLMETSIN[comps]"
-0100 %   can be used as inputs or outputs in the similar manner to normal
-0101 %   metabolites. This is a convenient way to, for example, allow excretion of
-0102 %   all metabolites to check whether it's the synthesis of some metabolite
-0103 %   that is limiting or whether it's the degradation of some byproduct. One
-0104 %   important difference is that only the upper bounds are used for these general
-0105 %   metabolites. That is, you can only say that uptake or excretion is
-0106 %   allowed, not that it is required. This is to avoid conflicts where the
-0107 %   constraints for the general metabolites overwrite those of the real
-0108 %   ones.
-0109 %
-0110 % Usage: taskStruct=parseTaskList(inputFile)
-0111 
-0112 if ~isfile(inputFile)
-0113     error('Task list %s cannot be found',string(inputFile));
-0114 end
-0115 
-0116 %Load the tasks file
-0117 convNumeric = false;
-0118 if strcmp(extractAfter(inputFile,strlength(inputFile) - 4), '.txt')
-0119     %load from tab delimited text file
-0120     fid = fopen(inputFile);
-0121     %Need to read numeric columns as strings, this is converted further
-0122     %down. If not, the titles would be lost.
-0123     convNumeric = true;
-0124     C = textscan(fid,'%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%*[^\n]', 'Delimiter', '\t');
-0125     fclose(fid);
-0126     raw = [C{:}];%unnest the cell array of cell arrays into a 2-dim cell array
-0127 else
-0128     [raw,flag]=loadSheet(loadWorkbook(inputFile), 'TASKS');
-0129     if flag~=0
-0130         EM=['Could not load sheet "TASKS" from ' inputFile];
-0131         dispEM(EM);
-0132     end
-0133 end
-0134 
-0135 %Remove all lines starting with "#" (or actually any character) and all
-0136 %empty columns
-0137 raw=cleanSheet(raw);
-0138 
-0139 %Captions
-0140 columns={'ID';'DESCRIPTION';'IN';'IN LB';'IN UB';'OUT';'OUT LB';'OUT UB';'EQU';'EQU LB';'EQU UB';'CHANGED RXN';'CHANGED LB';'CHANGED UB';'SHOULD FAIL';'PRINT FLUX';'COMMENTS'};
-0141 %Match the columns
-0142 [I, colI]=ismember(columns,raw(1,:));
-0143 
-0144 %If read from  a text file, the numbers will be strings - fix that
-0145 if convNumeric % in theory, this if should not be needed, the code should do nothing if all are already numeric. But it is kept as a safeguard.
-0146     numericColumns = [0 0 0 1 1 0 1 1 0 1 1 0 1 1 0 0 0] == 1;
-0147     cols = colI(numericColumns);
-0148     numeric = cellfun(@isnumeric,raw(:,cols));
-0149     %trick to avoid messing up the title row:
-0150     numeric(1,:) = 1;
-0151     for colind = 1:numel(cols)
-0152         col = cols(colind);
-0153         raw(~numeric(:,colind),col) = cellfun(@str2num, raw(~numeric(:,colind),col), 'UniformOutput', false);
-0154     end
-0155 end
-0156 
-0157 %Check that the ID field is present
-0158 if I(1)==0
-0159     EM='The TASKS sheet must have a column named ID';
-0160     dispEM(EM);
-0161 end
-0162 
-0163 %make sure numerical fields are converted from strings
-0164 
-0165 %Add default bounds where needed
-0166 for i=[4 5 7 8]
-0167     I=cellfun(@isempty,raw(:,colI(i)));
-0168     if i==5 || i==8
-0169         raw(I,colI(i))={1000};
-0170     else
-0171         raw(I,colI(i))={0};
-0172     end
-0173 end
-0174 
-0175 %Create an empty task structure
-0176 eTask.id='';
-0177 eTask.description='';
-0178 eTask.shouldFail=false;
-0179 eTask.printFluxes=false;
-0180 eTask.comments='';
-0181 eTask.inputs={};
-0182 eTask.LBin=[];
-0183 eTask.UBin=[];
-0184 eTask.outputs={};
-0185 eTask.LBout=[];
-0186 eTask.UBout=[];
-0187 eTask.equations={};
-0188 eTask.LBequ=[];
-0189 eTask.UBequ=[];
-0190 eTask.changed={};
-0191 eTask.LBrxn=[];
-0192 eTask.UBrxn=[];
-0193 
-0194 %Main loop
-0195 taskStruct=[];
-0196 task=eTask;
-0197 if isnumeric(raw{2,colI(1)})
-0198     task.id=num2str(raw{2,colI(1)});
-0199 else
-0200     task.id=raw{2,colI(1)};
-0201 end
-0202 task.description=raw{2,colI(2)};
-0203 if ~isempty(raw{2,colI(15)})
-0204     task.shouldFail=true;
-0205 end
-0206 if ~isempty(raw{2,colI(16)})
-0207     task.printFluxes=true;
-0208 end
-0209 if ~isempty(raw{2,colI(17)})
-0210     task.comments=raw{2,colI(17)};
-0211 end
-0212 
-0213 for i=2:size(raw,1)
-0214     %Set the inputs
-0215     if ischar(raw{i,colI(3)})
-0216         inputs=regexp(raw{i,colI(3)},';','split');
-0217         task.inputs=[task.inputs;inputs(:)];
-0218         task.LBin=[task.LBin;ones(numel(inputs),1)*raw{i,colI(4)}];
-0219         task.UBin=[task.UBin;ones(numel(inputs),1)*raw{i,colI(5)}];
-0220     end
-0221     %Set the outputs
-0222     if ischar(raw{i,colI(6)})
-0223         outputs=regexp(raw{i,colI(6)},';','split');
-0224         task.outputs=[task.outputs;outputs(:)];
-0225         task.LBout=[task.LBout;ones(numel(outputs),1)*raw{i,colI(7)}];
-0226         task.UBout=[task.UBout;ones(numel(outputs),1)*raw{i,colI(8)}];
-0227     end
-0228     %Add new rxns
-0229     if ischar(raw{i,colI(9)})
-0230         task.equations=[task.equations;raw{i,colI(9)}];
-0231         if ~isempty(raw{i,colI(10)})
-0232             task.LBequ=[task.LBequ;raw{i,colI(10)}];
-0233         else
-0234             if any(strfind(raw{i,colI(9)},'<=>'))
-0235                 task.LBequ=[task.LBequ;-1000];
-0236             else
-0237                 task.LBequ=[task.LBequ;0];
-0238             end
-0239         end
-0240         if ~isempty(raw{i,colI(11)})
-0241             task.UBequ=[task.UBequ;raw{i,colI(11)}];
-0242         else
-0243             task.UBequ=[task.UBequ;1000];
-0244         end
-0245     end
-0246     %Add changed bounds
-0247     if ischar(raw{i,colI(12)})
-0248         changed=regexp(raw{i,colI(12)},';','split');
-0249         task.changed=[task.changed;changed(:)];
-0250         task.LBrxn=[task.LBrxn;ones(numel(changed),1)*raw{i,colI(13)}];
-0251         task.UBrxn=[task.UBrxn;ones(numel(changed),1)*raw{i,colI(14)}];
-0252     end
-0253     
-0254     %Check if it should add more constraints
-0255     if i<size(raw,1)
-0256         if isempty(raw{i+1,colI(1)})
-0257             continue;
-0258         end
-0259     end
-0260     
-0261     taskStruct=[taskStruct;task];
-0262     task=eTask;
-0263     if i<size(raw,1)
-0264         if isnumeric(raw{i+1,colI(1)})
-0265             task.id=num2str(raw{i+1,colI(1)});
-0266         else
-0267             task.id=raw{i+1,colI(1)};
-0268         end
-0269         task.description=raw{i+1,colI(2)};
-0270         if ~isempty(raw{i+1,colI(15)})
-0271             task.shouldFail=true;
-0272         end
-0273         if ~isempty(raw{i+1,colI(16)})
-0274             task.printFluxes=true;
-0275         end
-0276         if ~isempty(raw{i+1,colI(17)})
-0277             task.comments=raw{i+1,colI(17)};
-0278         end
-0279     end
-0280 end
-0281 
-0282 %Should add more checks, such as unique IDs and missing headers
-0283 
-0284 end
+0005 % inputFile a task list in either Excel (*.xlsx, with a sheet named +0006 % TASKS with all relevant content) or tab-delimited +0007 % (*.txt) format. The file may contain the following +0008 % column headers (note, all rows starting with a +0009 % non-empty cell are removed. The first row after that +0010 % is considered the headers): +0011 % ID +0012 % the only required header. Each task must have a +0013 % unique id (string or numeric). Tasks can span multiple +0014 % rows, only the first row in each task should have +0015 % an id +0016 % DESCRIPTION +0017 % description of the task +0018 % IN +0019 % allowed input(s) for the task. Metabolite names +0020 % should be on the form +0021 % "model.metName[model.comps]". Several inputs +0022 % can be delimited by ";". If so, then the same +0023 % bounds are used for all inputs. If that is not +0024 % wanted, then use several rows for the task +0025 % IN LB +0026 % lower bound for the uptake of the metabolites in +0027 % the row (optional, default 0 which corresponds to a +0028 % minimal uptake of 0 units) +0029 % IN UB +0030 % upper bound for the uptake of the metabolites in +0031 % the row (optional, default 1000 which corresponds to a +0032 % maximal uptake of 1000 units) +0033 % OUT +0034 % allowed output(s) for the task (see IN) +0035 % OUT LB +0036 % lower bound for the production of the metabolites in +0037 % the row (optional, default 0 which corresponds to a +0038 % minimal production of 0 units) +0039 % OUT UB +0040 % upper bound for the production of the metabolites in +0041 % the row (optional, default 1000 which corresponds to a +0042 % maximal production of 1000 units) +0043 % EQU +0044 % equation to add. The equation should be on the form +0045 % "0.4 A + 2 B <=> (or =>) C" and the metabolites +0046 % should be on the form +0047 % "model.metName[model.comps]" (optional) +0048 % EQU LB +0049 % lower bound for the equation (optional, default -1000 +0050 % for reversible and 0 for irreversible) +0051 % EQU UB +0052 % upper bound for the equation (optional, default 1000) +0053 % CHANGED RXN +0054 % reaction ID for which to change the bounds for. +0055 % Several IDs can be delimited by ";". If so, +0056 % then the same bounds are used for all reactions. If +0057 % that is not wanted, then use several rows for the task +0058 % CHANGED LB +0059 % lower bound for the reaction +0060 % CHANGED UB +0061 % upper bound for the reaction +0062 % SHOULD FAIL +0063 % true if the correct behavior of the model is to +0064 % not have a feasible solution given the constraints +0065 % (optional, default false) +0066 % PRINT FLUX +0067 % true if the function should print the corresponding +0068 % flux distribution for a task. Can be useful for +0069 % testing (optional, default false) +0070 % +0071 % taskStruct array of structures with the following fields +0072 % id the id of the task +0073 % description the description of the task +0074 % shouldFail true if the task should fail +0075 % printFluxes true if the fluxes should be printed +0076 % comments string with comments +0077 % inputs cell array with input metabolites (in the form metName[comps]) +0078 % LBin array with lower bounds on inputs (default, 0) +0079 % UBin array with upper bounds on inputs (default, 1000) +0080 % outputs cell array with output metabolites (in the form metName[comps]) +0081 % LBout array with lower bounds on outputs (default, 0) +0082 % UBout array with upper bounds on outputs (default, 1000) +0083 % equations cell array with equations (with mets in the form metName[comps]) +0084 % LBequ array with lower bounds on equations (default, -1000 for +0085 % reversible and 0 for irreversible) +0086 % UBequ array with upper bounds on equations (default, 1000) +0087 % changed cell array with reactions to change bounds for +0088 % LBrxn array with lower bounds on changed reactions +0089 % UBrxn array with upper bounds on changed reactions +0090 % +0091 % This function is used for defining a set of tasks for a model to +0092 % perform. The tasks are defined by defining constraints on the model, +0093 % and if the problem is feasible, then the task is considered successful. +0094 % In general, each row can contain one constraint on uptakes, one +0095 % constraint on outputs, one new equation, and one change of reaction +0096 % bounds. If more bounds are needed to define the task, then several rows +0097 % can be used for each task. To perform the task use checkTasks or +0098 % fitTasks. +0099 % +0100 % NOTE: The general metabolites "ALLMETS" and "ALLMETSIN[comps]" +0101 % can be used as inputs or outputs in the similar manner to normal +0102 % metabolites. This is a convenient way to, for example, allow excretion of +0103 % all metabolites to check whether it's the synthesis of some metabolite +0104 % that is limiting or whether it's the degradation of some byproduct. One +0105 % important difference is that only the upper bounds are used for these general +0106 % metabolites. That is, you can only say that uptake or excretion is +0107 % allowed, not that it is required. This is to avoid conflicts where the +0108 % constraints for the general metabolites overwrite those of the real +0109 % ones. +0110 % +0111 % Usage: taskStruct=parseTaskList(inputFile) +0112 +0113 if ~isfile(inputFile) +0114 error('Task list %s cannot be found',string(inputFile)); +0115 end +0116 +0117 %Load the tasks file +0118 convNumeric = false; +0119 if strcmp(extractAfter(inputFile,strlength(inputFile) - 4), '.txt') +0120 %load from tab delimited text file +0121 fid = fopen(inputFile); +0122 %Need to read numeric columns as strings, this is converted further +0123 %down. If not, the titles would be lost. +0124 convNumeric = true; +0125 C = textscan(fid,'%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%q%*[^\n]', 'Delimiter', '\t'); +0126 fclose(fid); +0127 raw = [C{:}];%unnest the cell array of cell arrays into a 2-dim cell array +0128 else +0129 [raw,flag]=loadSheet(loadWorkbook(inputFile), 'TASKS'); +0130 if flag~=0 +0131 EM=['Could not load sheet "TASKS" from ' inputFile]; +0132 dispEM(EM); +0133 end +0134 end +0135 +0136 %Remove all lines starting with "#" (or actually any character) and all +0137 %empty columns +0138 raw=cleanSheet(raw); +0139 +0140 %Captions +0141 columns={'ID';'DESCRIPTION';'IN';'IN LB';'IN UB';'OUT';'OUT LB';'OUT UB';'EQU';'EQU LB';'EQU UB';'CHANGED RXN';'CHANGED LB';'CHANGED UB';'SHOULD FAIL';'PRINT FLUX';'COMMENTS'}; +0142 %Match the columns +0143 [I, colI]=ismember(columns,raw(1,:)); +0144 +0145 %If read from a text file, the numbers will be strings - fix that +0146 if convNumeric % in theory, this if should not be needed, the code should do nothing if all are already numeric. But it is kept as a safeguard. +0147 numericColumns = [0 0 0 1 1 0 1 1 0 1 1 0 1 1 0 0 0] == 1; +0148 cols = colI(numericColumns); +0149 numeric = cellfun(@isnumeric,raw(:,cols)); +0150 %trick to avoid messing up the title row: +0151 numeric(1,:) = 1; +0152 for colind = 1:numel(cols) +0153 col = cols(colind); +0154 raw(~numeric(:,colind),col) = cellfun(@str2num, raw(~numeric(:,colind),col), 'UniformOutput', false); +0155 end +0156 end +0157 +0158 %Check that the ID field is present +0159 if I(1)==0 +0160 EM='The TASKS sheet must have a column named ID'; +0161 dispEM(EM); +0162 end +0163 +0164 %make sure numerical fields are converted from strings +0165 +0166 %Add default bounds where needed +0167 for i=[4 5 7 8] +0168 I=cellfun(@isempty,raw(:,colI(i))); +0169 if i==5 || i==8 +0170 raw(I,colI(i))={1000}; +0171 else +0172 raw(I,colI(i))={0}; +0173 end +0174 end +0175 +0176 %Create an empty task structure +0177 eTask.id=''; +0178 eTask.description=''; +0179 eTask.shouldFail=false; +0180 eTask.printFluxes=false; +0181 eTask.comments=''; +0182 eTask.inputs={}; +0183 eTask.LBin=[]; +0184 eTask.UBin=[]; +0185 eTask.outputs={}; +0186 eTask.LBout=[]; +0187 eTask.UBout=[]; +0188 eTask.equations={}; +0189 eTask.LBequ=[]; +0190 eTask.UBequ=[]; +0191 eTask.changed={}; +0192 eTask.LBrxn=[]; +0193 eTask.UBrxn=[]; +0194 +0195 %Main loop +0196 taskStruct=[]; +0197 task=eTask; +0198 if isnumeric(raw{2,colI(1)}) +0199 task.id=num2str(raw{2,colI(1)}); +0200 else +0201 task.id=raw{2,colI(1)}; +0202 end +0203 task.description=raw{2,colI(2)}; +0204 if ~isempty(raw{2,colI(15)}) +0205 task.shouldFail=true; +0206 end +0207 if ~isempty(raw{2,colI(16)}) +0208 task.printFluxes=true; +0209 end +0210 if ~isempty(raw{2,colI(17)}) +0211 task.comments=raw{2,colI(17)}; +0212 end +0213 +0214 for i=2:size(raw,1) +0215 %Set the inputs +0216 if ischar(raw{i,colI(3)}) +0217 inputs=regexp(raw{i,colI(3)},';','split'); +0218 task.inputs=[task.inputs;inputs(:)]; +0219 task.LBin=[task.LBin;ones(numel(inputs),1)*raw{i,colI(4)}]; +0220 task.UBin=[task.UBin;ones(numel(inputs),1)*raw{i,colI(5)}]; +0221 end +0222 %Set the outputs +0223 if ischar(raw{i,colI(6)}) +0224 outputs=regexp(raw{i,colI(6)},';','split'); +0225 task.outputs=[task.outputs;outputs(:)]; +0226 task.LBout=[task.LBout;ones(numel(outputs),1)*raw{i,colI(7)}]; +0227 task.UBout=[task.UBout;ones(numel(outputs),1)*raw{i,colI(8)}]; +0228 end +0229 %Add new rxns +0230 if ischar(raw{i,colI(9)}) +0231 task.equations=[task.equations;raw{i,colI(9)}]; +0232 if ~isempty(raw{i,colI(10)}) +0233 task.LBequ=[task.LBequ;raw{i,colI(10)}]; +0234 else +0235 if any(strfind(raw{i,colI(9)},'<=>')) +0236 task.LBequ=[task.LBequ;-1000]; +0237 else +0238 task.LBequ=[task.LBequ;0]; +0239 end +0240 end +0241 if ~isempty(raw{i,colI(11)}) +0242 task.UBequ=[task.UBequ;raw{i,colI(11)}]; +0243 else +0244 task.UBequ=[task.UBequ;1000]; +0245 end +0246 end +0247 %Add changed bounds +0248 if ischar(raw{i,colI(12)}) +0249 changed=regexp(raw{i,colI(12)},';','split'); +0250 task.changed=[task.changed;changed(:)]; +0251 task.LBrxn=[task.LBrxn;ones(numel(changed),1)*raw{i,colI(13)}]; +0252 task.UBrxn=[task.UBrxn;ones(numel(changed),1)*raw{i,colI(14)}]; +0253 end +0254 +0255 %Check if it should add more constraints +0256 if i<size(raw,1) +0257 if isempty(raw{i+1,colI(1)}) +0258 continue; +0259 end +0260 end +0261 +0262 taskStruct=[taskStruct;task]; +0263 task=eTask; +0264 if i<size(raw,1) +0265 if isnumeric(raw{i+1,colI(1)}) +0266 task.id=num2str(raw{i+1,colI(1)}); +0267 else +0268 task.id=raw{i+1,colI(1)}; +0269 end +0270 task.description=raw{i+1,colI(2)}; +0271 if ~isempty(raw{i+1,colI(15)}) +0272 task.shouldFail=true; +0273 end +0274 if ~isempty(raw{i+1,colI(16)}) +0275 task.printFluxes=true; +0276 end +0277 if ~isempty(raw{i+1,colI(17)}) +0278 task.comments=raw{i+1,colI(17)}; +0279 end +0280 end +0281 end +0282 +0283 %Should add more checks, such as unique IDs and missing headers +0284 +0285 end
Generated by m2html © 2005
\ No newline at end of file From 7ced68b372cc9e11f8ccde9b2db539cb1b612b4f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 15 Jan 2026 22:43:21 +0100 Subject: [PATCH 6/9] fix: mergeModels and copyToComps ...From fields --- core/copyToComps.m | 40 +- core/mergeModels.m | 55 +- doc/core/copyToComps.html | 162 +-- doc/core/dispEM.html | 2 +- doc/core/fillGaps.html | 6 +- doc/core/fitTasks.html | 4 +- doc/core/getModelFromHomology.html | 4 +- doc/core/mergeModels.html | 1239 +++++++++-------- doc/core/standardizeGrRules.html | 2 +- doc/index.html | 113 +- doc/testing/unit_tests/importExportTests.html | 2 +- doc/tutorial/index.html | 4 +- doc/utils/emptyOrLogicalScalar.html | 53 + doc/utils/emptyOrTextOrCellOfText.html | 90 ++ doc/utils/emptyOrTextScalar.html | 59 + doc/utils/index.html | 29 + testing/unit_tests/importExportTests.m | 2 +- utils/emptyOrLogicalScalar.m | 6 + utils/emptyOrTextOrCellOfText.m | 43 + utils/emptyOrTextScalar.m | 12 + 20 files changed, 1123 insertions(+), 804 deletions(-) create mode 100644 doc/utils/emptyOrLogicalScalar.html create mode 100644 doc/utils/emptyOrTextOrCellOfText.html create mode 100644 doc/utils/emptyOrTextScalar.html create mode 100644 doc/utils/index.html create mode 100644 utils/emptyOrLogicalScalar.m create mode 100644 utils/emptyOrTextOrCellOfText.m create mode 100644 utils/emptyOrTextScalar.m diff --git a/core/copyToComps.m b/core/copyToComps.m index 37dd4ccb..042b4c0b 100755 --- a/core/copyToComps.m +++ b/core/copyToComps.m @@ -25,24 +25,26 @@ % % Usage: model=copyToComps(model,toComps,rxns,deleteOriginal,compNames,compOutside) -if nargin<3 - rxns=model.rxns; -elseif ~islogical(rxns) && ~isnumeric(rxns) - rxns=convertCharArray(rxns); +arguments + model (1,1) struct + toComps {emptyOrTextOrCellOfText} + rxns = model.rxns + deleteOriginal {emptyOrLogicalScalar} = false + compNames {emptyOrTextOrCellOfText} = toComps + compOutside {emptyOrTextOrCellOfText} = ''; end -if nargin<4 - deleteOriginal=false; + +if nargin >= 3 && ~islogical(rxns) && ~isnumeric(rxns) + rxns = convertCharArray(rxns); end -if nargin<5 - compNames=toComps; -else +if nargin >= 5 compNames=convertCharArray(compNames); end -if nargin<6 - compOutside=cell(numel(toComps),1); - compOutside(:)={''}; -else +if nargin >= 6 compOutside=convertCharArray(compOutside); + if length(compOutside) ~= length(compNames) + error('compOutside and compNames should be of equal size.'); + end end originalID=model.id; @@ -90,17 +92,7 @@ end %Merge the models - model=mergeModels({model;modelToAdd},'metNames'); -end - -if all(strcmp(model.rxnFrom,originalID)) - model = rmfield(model,'rxnFrom'); -end -if all(strcmp(model.metFrom,originalID)) - model = rmfield(model,'metFrom'); -end -if all(strcmp(model.geneFrom,originalID)) - model = rmfield(model,'geneFrom'); + model=mergeModels({model;modelToAdd},'metNames',[],true); end if deleteOriginal==true diff --git a/core/mergeModels.m b/core/mergeModels.m index 47ad5007..cf45dd7a 100755 --- a/core/mergeModels.m +++ b/core/mergeModels.m @@ -1,4 +1,4 @@ -function model=mergeModels(models,metParam,supressWarnings) +function model=mergeModels(models,metParam,supressWarnings,copyToComps) % mergeModels % Merges models into one model structure. Reactions are added without any % checks, so duplicate reactions might appear. Metabolites are matched by @@ -8,9 +8,11 @@ % Input: % models a cell array with model structures % metParam string metabolite name ('metNames') or ID ('mets') are -% used for matching (default, metNames) +% 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 @@ -22,38 +24,49 @@ % % Usage: model=mergeModels(models) +arguments + models; + metParam {emptyOrTextScalar} = "metNames" + supressWarnings {emptyOrLogicalScalar} = false + copyToComps {emptyOrLogicalScalar} = false +end + +metParam = char(metParam); + %Just return the model if numel(models)<=1 model=models{1}; return; end -if nargin<2 - metParam='metNames'; -else - metParam=char(metParam); -end - -if nargin<3 - supressWarnings=false; -end - hasMetFrom = cellfun(@(s) isfield(s,'metFrom'), models); hasGeneFrom = cellfun(@(s) isfield(s,'geneFrom'), models); hasRxnFrom = cellfun(@(s) isfield(s,'rxnFrom'), models); for i = 1:numel(models) - if ~any(hasMetFrom) + if copyToComps + if hasMetFrom(1) + models{2}.metFrom = repmat({''},numel(models{i}.mets),1); + end + elseif ~any(hasMetFrom) models{i}.metFrom = repmat({models{i}.id},numel(models{i}.mets),1); elseif ~hasMetFrom(i) models{i}.metFrom = repmat({''},numel(models{i}.mets),1); end - if ~any(hasRxnFrom) + if copyToComps + if hasRxnFrom(1) + models{2}.rxnFrom = repmat({''},numel(models{i}.rxns),1); + end + elseif ~any(hasRxnFrom) models{i}.rxnFrom = repmat({models{i}.id},numel(models{i}.rxns),1); elseif ~hasRxnFrom(i) models{i}.rxnFrom = repmat({''},numel(models{i}.rxns),1); end - if ~any(hasGeneFrom) && any(cellfun(@(s) isfield(s,'genes'), models)) + if copyToComps + if hasGeneFrom(1) + models{2}.geneFrom = repmat({''},numel(models{i}.genes),1); + end + elseif ~any(hasGeneFrom) && any(cellfun(@(s) isfield(s,'genes'), models)) models{i}.geneFrom = repmat({models{i}.id},numel(models{i}.genes),1); elseif ~hasGeneFrom(i) models{i}.geneFrom = repmat({''},numel(models{i}.genes),1); @@ -95,7 +108,9 @@ end %Add all static stuff - model.rxnFrom = [model.rxnFrom; models{i}.rxnFrom]; + if any(hasRxnFrom) || (~copyToComps && ~any(hasRxnFrom)) + model.rxnFrom = [model.rxnFrom; models{i}.rxnFrom]; + end model.rxns = [model.rxns; models{i}.rxns]; model.rxnNames = [model.rxnNames; models{i}.rxnNames]; model.lb = [model.lb; models{i}.lb]; @@ -302,7 +317,9 @@ end %Add static info on the metabolites - model.metFrom = [model.metFrom; models{i}.metFrom(metsToAdd)]; + if any(hasMetFrom) + model.metFrom = [model.metFrom; models{i}.metFrom(metsToAdd)]; + end model.mets = [model.mets; models{i}.mets(metsToAdd)]; model.metNames = [model.metNames; models{i}.metNames(metsToAdd)]; model.b = [model.b; zeros(numel(metsToAdd),size(model.b,2))]; @@ -498,7 +515,9 @@ model.rxnGeneMat = [sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat]; emptyGene = repmat({''},numel(model.rxns),1); model.grRules = [emptyGene;models{i}.grRules]; - model.geneFrom = models{i}.geneFrom; + if any(hasGeneFrom) + model.geneFrom = models{i}.geneFrom; + end if isfield(models{i},'geneShortNames') model.geneShortNames=models{i}.geneShortNames; diff --git a/doc/core/copyToComps.html b/doc/core/copyToComps.html index 20c77753..33b17a6e 100644 --- a/doc/core/copyToComps.html +++ b/doc/core/copyToComps.html @@ -57,10 +57,10 @@

DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

This function calls: +
  • convertCharArray convertCharArray
  • getIndexes getIndexes
  • mergeModels mergeModels
  • removeReactions removeReactions
    • This function is called by:
      -
    +
  • mergeModels mergeModels
    • @@ -93,91 +93,83 @@

    SOURCE CODE ^% 0026 % Usage: model=copyToComps(model,toComps,rxns,deleteOriginal,compNames,compOutside) 0027 -0028 if nargin<3 -0029 rxns=model.rxns; -0030 elseif ~islogical(rxns) && ~isnumeric(rxns) -0031 rxns=convertCharArray(rxns); -0032 end -0033 if nargin<4 -0034 deleteOriginal=false; +0028 arguments +0029 model (1,1) struct +0030 toComps {emptyOrTextOrCellOfText} +0031 rxns = model.rxns +0032 deleteOriginal {emptyOrLogicalScalar} = false +0033 compNames {emptyOrTextOrCellOfText} = toComps +0034 compOutside {emptyOrTextOrCellOfText} = ''; 0035 end -0036 if nargin<5 -0037 compNames=toComps; -0038 else -0039 compNames=convertCharArray(compNames); -0040 end -0041 if nargin<6 -0042 compOutside=cell(numel(toComps),1); -0043 compOutside(:)={''}; -0044 else -0045 compOutside=convertCharArray(compOutside); -0046 end -0047 -0048 originalID=model.id; -0049 originalName=model.name; -0050 -0051 rxns=getIndexes(model,rxns,'rxns'); +0036 +0037 if nargin >= 3 && ~islogical(rxns) && ~isnumeric(rxns) +0038 rxns = convertCharArray(rxns); +0039 end +0040 if nargin >= 5 +0041 compNames=convertCharArray(compNames); +0042 end +0043 if nargin >= 6 +0044 compOutside=convertCharArray(compOutside); +0045 if length(compOutside) ~= length(compNames) +0046 error('compOutside and compNames should be of equal size.'); +0047 end +0048 end +0049 +0050 originalID=model.id; +0051 originalName=model.name; 0052 -0053 for i=1:numel(toComps) -0054 %Check if the compartment exists, otherwise add it -0055 [I,J]=ismember(toComps(i),model.comps); -0056 if I==false -0057 model.comps=[model.comps;toComps(i)]; -0058 model.compNames=[model.compNames;compNames(i)]; -0059 if isfield(model,'compOutside') -0060 model.compOutside=[model.compOutside;compOutside(i)]; -0061 end -0062 if isfield(model,'compMiriams') -0063 model.compMiriams=[model.compMiriams;cell(1,1)]; -0064 end -0065 J=numel(model.comps); -0066 end -0067 %Copy the reactions by making a model structure with only them, then -0068 %change the localization, and finally merge with the original model -0069 modelToAdd=model; -0070 modelToAdd=removeReactions(modelToAdd,setdiff(1:numel(model.rxns),rxns),true,true); -0071 modelToAdd.rxns=strcat(modelToAdd.rxns,'_',toComps(i)); -0072 modelToAdd.mets=strcat(modelToAdd.mets,'_',toComps(i)); -0073 modelToAdd.comps=modelToAdd.comps(J); -0074 modelToAdd.compNames=modelToAdd.compNames(J); -0075 if isfield(modelToAdd,'compOutside') -0076 modelToAdd.compOutside=modelToAdd.compOutside(J); -0077 end -0078 if isfield(modelToAdd,'compMiriams') -0079 modelToAdd.compMiriams=modelToAdd.compMiriams(J); -0080 end -0081 modelToAdd.metComps=ones(numel(modelToAdd.mets),1); -0082 if isfield(modelToAdd,'metFrom') -0083 modelToAdd = rmfield(modelToAdd,'metFrom'); -0084 end -0085 if isfield(modelToAdd,'rxnFrom') -0086 modelToAdd = rmfield(modelToAdd,'rxnFrom'); -0087 end -0088 if isfield(modelToAdd,'geneFrom') -0089 modelToAdd = rmfield(modelToAdd,'geneFrom'); -0090 end -0091 -0092 %Merge the models -0093 model=mergeModels({model;modelToAdd},'metNames'); -0094 end -0095 -0096 if all(strcmp(model.rxnFrom,originalID)) -0097 model = rmfield(model,'rxnFrom'); -0098 end -0099 if all(strcmp(model.metFrom,originalID)) -0100 model = rmfield(model,'metFrom'); -0101 end -0102 if all(strcmp(model.geneFrom,originalID)) -0103 model = rmfield(model,'geneFrom'); -0104 end -0105 -0106 if deleteOriginal==true -0107 model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments -0108 end -0109 -0110 model.id=originalID; -0111 model.name=originalName; -0112 end +0053 rxns=getIndexes(model,rxns,'rxns'); +0054 +0055 for i=1:numel(toComps) +0056 %Check if the compartment exists, otherwise add it +0057 [I,J]=ismember(toComps(i),model.comps); +0058 if I==false +0059 model.comps=[model.comps;toComps(i)]; +0060 model.compNames=[model.compNames;compNames(i)]; +0061 if isfield(model,'compOutside') +0062 model.compOutside=[model.compOutside;compOutside(i)]; +0063 end +0064 if isfield(model,'compMiriams') +0065 model.compMiriams=[model.compMiriams;cell(1,1)]; +0066 end +0067 J=numel(model.comps); +0068 end +0069 %Copy the reactions by making a model structure with only them, then +0070 %change the localization, and finally merge with the original model +0071 modelToAdd=model; +0072 modelToAdd=removeReactions(modelToAdd,setdiff(1:numel(model.rxns),rxns),true,true); +0073 modelToAdd.rxns=strcat(modelToAdd.rxns,'_',toComps(i)); +0074 modelToAdd.mets=strcat(modelToAdd.mets,'_',toComps(i)); +0075 modelToAdd.comps=modelToAdd.comps(J); +0076 modelToAdd.compNames=modelToAdd.compNames(J); +0077 if isfield(modelToAdd,'compOutside') +0078 modelToAdd.compOutside=modelToAdd.compOutside(J); +0079 end +0080 if isfield(modelToAdd,'compMiriams') +0081 modelToAdd.compMiriams=modelToAdd.compMiriams(J); +0082 end +0083 modelToAdd.metComps=ones(numel(modelToAdd.mets),1); +0084 if isfield(modelToAdd,'metFrom') +0085 modelToAdd = rmfield(modelToAdd,'metFrom'); +0086 end +0087 if isfield(modelToAdd,'rxnFrom') +0088 modelToAdd = rmfield(modelToAdd,'rxnFrom'); +0089 end +0090 if isfield(modelToAdd,'geneFrom') +0091 modelToAdd = rmfield(modelToAdd,'geneFrom'); +0092 end +0093 +0094 %Merge the models +0095 model=mergeModels({model;modelToAdd},'metNames',[],true); +0096 end +0097 +0098 if deleteOriginal==true +0099 model=removeReactions(model,rxns,true,true,true); %Also delete unused compartments +0100 end +0101 +0102 model.id=originalID; +0103 model.name=originalName; +0104 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/core/dispEM.html b/doc/core/dispEM.html index cd276912..745229f2 100644 --- a/doc/core/dispEM.html +++ b/doc/core/dispEM.html @@ -49,7 +49,7 @@

    CROSS-REFERENCE INFORMATION ^ <li><a href=convertCharArray convertCharArray This function is called by: +
  • addGenesRaven addGenesRaven
  • addMets addMets
  • addRxns addRxns
  • addRxnsGenesMets addRxnsGenesMets
  • addTransport addTransport
  • analyzeSampling analyzeSampling
  • buildEquation buildEquation
  • changeRxns changeRxns
  • checkModelStruct checkModelStruct
  • checkRxn checkRxn
  • checkTasks checkTasks
  • compareMultipleModels compareMultipleModels
  • compareRxnsGenesMetsComps compareRxnsGenesMetsComps
  • constructS constructS
  • consumeSomething consumeSomething
  • contractModel contractModel
  • fillGaps fillGaps
  • findGeneDeletions findGeneDeletions
  • fitParameters fitParameters
  • fitTasks fitTasks
  • getElementalBalance getElementalBalance
  • getEssentialRxns getEssentialRxns
  • getExpressionStructure getExpressionStructure
  • getFluxZ getFluxZ
  • getMetsInComp getMetsInComp
  • getMinNrFluxes getMinNrFluxes
  • getModelFromHomology getModelFromHomology
  • getRxnsInComp getRxnsInComp
  • guessComposition guessComposition
  • makeSomething makeSomething
  • mapCompartments mapCompartments
  • mergeCompartments mergeCompartments
  • mergeModels mergeModels
  • parseTaskList parseTaskList
  • predictLocalization predictLocalization
  • printFluxes printFluxes
  • randomSampling randomSampling
  • removeBadRxns removeBadRxns
  • reporterMetabolites reporterMetabolites
  • setParam setParam
  • simplifyModel simplifyModel
  • sortModel sortModel
    • diff --git a/doc/core/fillGaps.html b/doc/core/fillGaps.html index 5ff4d82e..ba68613a 100644 --- a/doc/core/fillGaps.html +++ b/doc/core/fillGaps.html @@ -94,7 +94,7 @@

    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

    This function calls: +
  • dispEM dispEM
  • getExchangeRxns getExchangeRxns
  • getMinNrFluxes getMinNrFluxes
  • haveFlux haveFlux
  • mergeModels mergeModels
  • removeReactions removeReactions
  • simplifyModel simplifyModel
    • This function is called by: @@ -235,7 +235,7 @@

    SOURCE CODE ^%First merge all models into one big one -0133 allModels=mergeModels([{model};models],'metNames',true); +0133 allModels=mergeModels([{model};models],'metNames',true); 0134 0135 %Add that net production is ok 0136 if allowNetProduction==true @@ -316,7 +316,7 @@

    SOURCE CODE ^removeReactions(allModels,I,true,true,true); 0213 -0214 newModel=mergeModels({model;addedModel},'metNames',true); +0214 newModel=mergeModels({model;addedModel},'metNames',true); 0215 addedRxns=setdiff(newModel.rxns,model.rxns); 0216 newModel=rmfield(newModel,'rxnScores'); 0217 end diff --git a/doc/core/fitTasks.html b/doc/core/fitTasks.html index fa51abdd..3a45616c 100644 --- a/doc/core/fitTasks.html +++ b/doc/core/fitTasks.html @@ -68,7 +68,7 @@

    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

    This function calls: +
  • addMets addMets
  • addRxns addRxns
  • dispEM dispEM
  • fillGaps fillGaps
  • mergeModels mergeModels
  • parseTaskList parseTaskList
  • printFluxes printFluxes
  • removeReactions removeReactions
  • setParam setParam
    • This function is called by:
    @@ -362,7 +362,7 @@

    SOURCE CODE ^%Add the reactions to the base model. It is not correct 0284 %to use newModel directly, as it may contain 0285 %reactions/constraints that are specific to this task -0286 model=mergeModels({model,removeReactions(newModel,setdiff(newModel.rxns,newRxns),true,true)},'metNames',true); +0286 model=mergeModels({model,removeReactions(newModel,setdiff(newModel.rxns,newRxns),true,true)},'metNames',true); 0287 0288 %Keep track of the added reactions 0289 addedRxns(ismember(refModel.rxns,newRxns),i)=true; diff --git a/doc/core/getModelFromHomology.html b/doc/core/getModelFromHomology.html index e8f6d374..5291424f 100644 --- a/doc/core/getModelFromHomology.html +++ b/doc/core/getModelFromHomology.html @@ -96,7 +96,7 @@

    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

    This function calls: +
  • convertCharArray convertCharArray
  • deleteUnusedGenes deleteUnusedGenes
  • dispEM dispEM
  • mergeModels mergeModels
  • removeGenes removeGenes
  • removeReactions removeReactions
  • standardizeGrRules standardizeGrRules
    • This function is called by:
    @@ -631,7 +631,7 @@

    SOURCE CODE ^%Now merge the models. All information should be correct except for 'or' 0526 %complexes -0527 draftModel=mergeModels(models,'metNames'); +0527 draftModel=mergeModels(models,'metNames'); 0528 0529 %Remove unnecessary OLD_ genes, that were added with OR relationships 0530 regexStr=['OLD_(', strjoin(modelNames(:),'|'),')_(\S^\))+']; diff --git a/doc/core/mergeModels.html b/doc/core/mergeModels.html index 530ecf62..bdae72a3 100644 --- a/doc/core/mergeModels.html +++ b/doc/core/mergeModels.html @@ -24,7 +24,7 @@

    PURPOSE ^mergeModels

    SYNOPSIS ^

    -
    function model=mergeModels(models,metParam,supressWarnings)
    +
    function model=mergeModels(models,metParam,supressWarnings,copyToComps)

    DESCRIPTION ^

     mergeModels
@@ -36,9 +36,11 @@ 
    DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
    
 This function calls:
 
+
  • copyToComps	copyToComps
  • dispEM	dispEM
  • standardizeGrRules	standardizeGrRules
    • 
 This function is called by:
 
@@ -63,7 +65,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 
 
 <h2><a name=SOURCE CODE ^
    
-
    0001 function model=mergeModels(models,metParam,supressWarnings)
+
    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
@@ -73,617 +75,636 @@ 
    SOURCE CODE ^% Input:
 0009 %   models          a cell array with model structures
 0010 %   metParam        string metabolite name ('metNames') or ID ('mets') are
-0011 %                   used for matching (default, metNames)
+0011 %                   used for matching (optional, default 'metNames')
 0012 %   supressWarnings logical whether warnings should be supressed (optional,
 0013 %                   default false)
-0014 %
-0015 % Output:
-0016 %   model           a model structure with the merged model. Follows the
-0017 %                   structure of normal models but also has 'rxnFrom/
-0018 %                   metFrom/geneFrom' fields to indicate from which model
-0019 %                   each reaction/metabolite/gene was taken. If the model
-0020 %                   already has 'rxnFrom/metFrom/geneFrom' fields, then
-0021 %                   these fields are not modified.
-0022 %
-0023 % Usage: model=mergeModels(models)
-0024 
-0025 %Just return the model
-0026 if numel(models)<=1
-0027     model=models{1};
-0028     return;
-0029 end
-0030 
-0031 if nargin<2
-0032     metParam='metNames';
-0033 else
-0034     metParam=char(metParam);
-0035 end
-0036 
-0037 if nargin<3
-0038     supressWarnings=false;
-0039 end
-0040 
-0041 hasMetFrom  = cellfun(@(s) isfield(s,'metFrom'), models);
-0042 hasGeneFrom = cellfun(@(s) isfield(s,'geneFrom'), models);
-0043 hasRxnFrom  = cellfun(@(s) isfield(s,'rxnFrom'), models);
-0044 
-0045 for i = 1:numel(models)
-0046     if ~any(hasMetFrom)
-0047         models{i}.metFrom = repmat({models{i}.id},numel(models{i}.mets),1);
-0048     elseif ~hasMetFrom(i)
-0049         models{i}.metFrom = repmat({''},numel(models{i}.mets),1);
-0050     end
-0051     if ~any(hasRxnFrom)
-0052         models{i}.rxnFrom = repmat({models{i}.id},numel(models{i}.rxns),1);
-0053     elseif ~hasRxnFrom(i)
-0054         models{i}.rxnFrom = repmat({''},numel(models{i}.rxns),1);
+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 ~any(hasGeneFrom) && any(cellfun(@(s) isfield(s,'genes'), models))
-0057         models{i}.geneFrom = repmat({models{i}.id},numel(models{i}.genes),1);
-0058     elseif ~hasGeneFrom(i)
-0059         models{i}.geneFrom = repmat({''},numel(models{i}.genes),1);
-0060     end
-0061 end
-0062 
-0063 %Add new functionality in the order specified in models
-0064 model=models{1};
-0065 model.id='MERGED';
-0066 model.name='';
-0067 
-0068 if isfield(model,'equations')
-0069     model=rmfield(model,'equations');
-0070 end
-0071 
-0072 for i=2:numel(models)
-0073     %Add the model id to the rxn id id it already exists in the model (id
-0074     %have to be unique) This is because it makes a '[]' string if no new
-0075     %reactions
-0076     if ~isempty(models{i}.rxns)
-0077         I=ismember(models{i}.rxns,model.rxns);
-0078         models{i}.rxns(I)=strcat(models{i}.rxns(I),['_' models{i}.id]);
-0079     end
-0080     
-0081     %Make sure that there are no conflicting reaction ids
-0082     [~, ~, conflicting]=intersect(model.rxns,models{i}.rxns);
-0083     
-0084     if ~isempty(conflicting)
-0085         printString=cell(numel(conflicting),1);
-0086         for j=1:numel(conflicting)
-0087             printString{j}=['Old: ' models{i}.rxns{conflicting(j)} ' New: ' models{i}.rxns{conflicting(j)} '_' models{i}.id];
-0088             models{i}.rxns{conflicting(j)}=[models{i}.rxns{conflicting(j)} '_' models{i}.id];
-0089         end
-0090         if supressWarnings==false
-0091             EM=['The following reaction IDs in ' models{i}.id ' are already present in the model and were renamed:'];
-0092             dispEM(EM,false,printString);
-0093             fprintf('\n');
-0094         end
-0095     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=2:numel(models)
+0086     %Add the model id to the rxn id id it already exists in the model (id
+0087     %have to be unique) This is because it makes a '[]' string if no new
+0088     %reactions
+0089     if ~isempty(models{i}.rxns)
+0090         I=ismember(models{i}.rxns,model.rxns);
+0091         models{i}.rxns(I)=strcat(models{i}.rxns(I),['_' models{i}.id]);
+0092     end
+0093     
+0094     %Make sure that there are no conflicting reaction ids
+0095     [~, ~, conflicting]=intersect(model.rxns,models{i}.rxns);
 0096     
-0097     %Add all static stuff
-0098     model.rxnFrom  = [model.rxnFrom;  models{i}.rxnFrom];
-0099     model.rxns     = [model.rxns;     models{i}.rxns];
-0100     model.rxnNames = [model.rxnNames; models{i}.rxnNames];
-0101     model.lb       = [model.lb;       models{i}.lb];
-0102     model.ub       = [model.ub;       models{i}.ub];
-0103     model.c        = [model.c;        models{i}.c];
-0104     model.rev      = [model.rev;      models{i}.rev];
-0105     
-0106     if isfield(models{i},'subSystems')
-0107         if isfield(model,'subSystems')
-0108             model.subSystems=[model.subSystems;models{i}.subSystems];
-0109         else
-0110             emptySubSystem=cell(numel(model.rxns)-numel(models{i}.rxns),1);
-0111             emptySubSystem(:)={''};
-0112             emptySubSystem=cellfun(@(x) cell(0,0),emptySubSystem,'UniformOutput',false);
-0113             model.subSystems=[emptySubSystem;models{i}.subSystems];
-0114         end
-0115     else
-0116         if isfield(model,'subSystems')
-0117             emptySubSystem=cell(numel(models{i}.rxns),1);
-0118             emptySubSystem(:)={''};
-0119             emptySubSystem=cellfun(@(x) cell(0,0),emptySubSystem,'UniformOutput',false);
-0120             model.subSystems=[model.subSystems;emptySubSystem];
-0121         end
-0122     end
-0123     
-0124     if isfield(models{i},'eccodes')
-0125         if isfield(model,'eccodes')
-0126             model.eccodes=[model.eccodes;models{i}.eccodes];
-0127         else
-0128             emptyEC=cell(numel(model.rxns)-numel(models{i}.rxns),1);
-0129             emptyEC(:)={''};
-0130             model.eccodes=[emptyEC;models{i}.eccodes];
-0131         end
-0132     else
-0133         if isfield(model,'eccodes')
-0134             emptyEC=cell(numel(models{i}.rxns),1);
-0135             emptyEC(:)={''};
-0136             model.eccodes=[model.eccodes;emptyEC];
-0137         end
-0138     end
-0139     
-0140     if isfield(models{i},'rxnMiriams')
-0141         if isfield(model,'rxnMiriams')
-0142             model.rxnMiriams=[model.rxnMiriams;models{i}.rxnMiriams];
-0143         else
-0144             model.rxnMiriams=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnMiriams];
-0145         end
-0146     else
-0147         if isfield(model,'rxnMiriams')
-0148             model.rxnMiriams=[model.rxnMiriams;cell(numel(models{i}.rxns),1)];
-0149         end
-0150     end
-0151     
-0152     if isfield(models{i},'rxnNotes')
-0153         if isfield(model,'rxnNotes')
-0154             model.rxnNotes=[model.rxnNotes;models{i}.rxnNotes];
-0155         else
-0156             emptyNotes=cell(numel(model.rxns)-numel(models{i}.rxns),1);
-0157             emptyNotes(:)={''};
-0158             model.rxnNotes=[emptyNotes;models{i}.rxnNotes];
-0159         end
-0160     else
-0161         if isfield(model,'rxnNotes')
-0162             emptyNotes=cell(numel(models{i}.rxns),1);
-0163             emptyNotes(:)={''};
-0164             model.rxnNotes=[model.rxnNotes;emptyNotes];
-0165         end
-0166     end
-0167     
-0168     if isfield(models{i},'rxnReferences')
-0169         if isfield(model,'rxnReferences')
-0170             model.rxnReferences=[model.rxnReferences;models{i}.rxnReferences];
-0171         else
-0172             emptyReferences=cell(numel(model.rxns)-numel(models{i}.rxns),1);
-0173             emptyReferences(:)={''};
-0174             model.rxnReferences=[emptyReferences;models{i}.rxnReferences];
-0175         end
-0176     else
-0177         if isfield(model,'rxnReferences')
-0178             emptyReferences=cell(numel(models{i}.rxns),1);
-0179             emptyReferences(:)={''};
-0180             model.rxnReferences=[model.rxnReferences;emptyReferences];
-0181         end
-0182     end
-0183     
-0184     if isfield(models{i},'rxnConfidenceScores')
-0185         if isfield(model,'rxnConfidenceScores')
-0186             model.rxnConfidenceScores=[model.rxnConfidenceScores;models{i}.rxnConfidenceScores];
-0187         else
-0188             model.rxnConfidenceScores=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnConfidenceScores];
-0189         end
-0190     else
-0191         if isfield(model,'rxnConfidenceScores')
-0192             model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(models{i}.rxns),1)];
-0193         end
-0194     end
-0195 
-0196     if isfield(models{i},'rxnDeltaG')
-0197         if isfield(model,'rxnDeltaG')
-0198             model.rxnDeltaG=[model.rxnDeltaG;models{i}.rxnDeltaG];
-0199         else
-0200             model.rxnDeltaG=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnDeltaG];
-0201         end
-0202     else
-0203         if isfield(model,'rxnDeltaG')
-0204             model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(models{i}.rxns),1)];
-0205         end
-0206     end
-0207     
-0208     if isfield(models{i},'rxnComps')
-0209         if isfield(model,'rxnComps')
-0210             model.rxnComps=[model.rxnComps;models{i}.rxnComps];
-0211         else
-0212             model.rxnComps=[ones(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnComps];
-0213             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');
-0214         end
-0215     else
-0216         if isfield(model,'rxnComps')
-0217             model.rxnComps=[model.rxnComps;ones(numel(models{i}.rxns),1)];
-0218             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');
-0219         end
-0220     end
-0221     
-0222     if isfield(models{i},'rxnScores')
-0223         if isfield(model,'rxnScores')
-0224             model.rxnScores=[model.rxnScores;models{i}.rxnScores];
-0225         else
-0226             emptyRS=zeros(numel(model.rxns)-numel(models{i}.rxns),1);
-0227             model.rxnScores=[emptyRS;models{i}.rxnScores];
-0228         end
-0229     else
-0230         if isfield(model,'rxnScores')
-0231             emptyRS=zeros(numel(models{i}.rxns),1);
-0232             model.rxnScores=[model.rxnScores;emptyRS];
-0233         end
-0234     end
-0235     
-0236     if isfield(models{i},'pwys')
-0237         if isfield(model,'pwys')
-0238             model.pwys=[model.pwys;models{i}.pwys];
-0239         else
-0240             model.pwys=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.pwys];
-0241         end
-0242     else
-0243         if isfield(model,'pwys')
-0244             model.pwys=[model.pwys;cell(numel(models{i}.rxns),1)];
-0245         end
-0246     end
-0247 
-0248     if strcmpi(metParam,'metNames')
-0249     %Get the new metabolites from matching the models. Metabolites are said
-0250     %to be the same if they share name and compartment id. This means that
-0251     %metabolite IDs are not taken into account.
-0252         
-0253         oldMetComps=model.comps(model.metComps);
-0254         oldMets=strcat(model.metNames,'[',oldMetComps,']');
-0255         %This is because it makes a '[]' string if no new metabolites
-0256         if ~isempty(models{i}.metNames)
-0257             newMetComps=models{i}.comps(models{i}.metComps);
-0258             newMets=strcat(models{i}.metNames,'[',newMetComps,']');
-0259         else
-0260             newMets={};
-0261             newMetComps={};
-0262         end
-0263         tf=ismember(newMets,oldMets);
-0264         metsToAdd=find(~tf);
-0265 
-0266     end
-0267 
-0268     if strcmpi(metParam,'mets')
-0269     %Get the new metabolites from matching the models. Metabolites are matched by metabolite ID (model.mets).
-0270 
-0271         oldMetComps=model.comps(model.metComps);
-0272         oldMets=model.mets;
-0273     
-0274         if ~isempty(models{i}.mets)
-0275             newMetComps=models{i}.comps(models{i}.metComps);
-0276             newMets=models{i}.mets;
-0277         else
-0278             newMets={};
-0279             newMetComps={};
-0280         end
-0281         tf=ismember(newMets,oldMets);
-0282         metsToAdd=find(~tf);
-0283 
-0284     end
-0285     
-0286     %First add the new metabolites Make sure that there are no conflicting
-0287     %metabolite ids
-0288     conflicting=ismember(models{i}.mets(metsToAdd),model.mets);
-0289     
-0290     conflicting=find(conflicting);
-0291     
-0292     if ~isempty(conflicting)
-0293         printString=cell(numel(conflicting),1);
-0294         for j=1:numel(conflicting)
-0295             printString{j}=['Old: ' models{i}.mets{metsToAdd(conflicting(j))} ' New: ' models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
-0296             models{i}.mets{metsToAdd(conflicting(j))}=[models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
-0297         end
-0298         if supressWarnings==false
-0299             EM=['The following metabolite IDs in ' models{i}.id ' are already present in the model and were renamed:'];
-0300             dispEM(EM,false,printString);
-0301         end
-0302     end
-0303     
-0304     %Add static info on the metabolites
-0305     model.metFrom  = [model.metFrom;  models{i}.metFrom(metsToAdd)];
-0306     model.mets     = [model.mets;     models{i}.mets(metsToAdd)];
-0307     model.metNames = [model.metNames; models{i}.metNames(metsToAdd)];
-0308     model.b        = [model.b;        zeros(numel(metsToAdd),size(model.b,2))];
-0309     
-0310     if isfield(model,'unconstrained')
-0311         if isfield(models{i},'unconstrained')
-0312             model.unconstrained=[model.unconstrained;models{i}.unconstrained(metsToAdd)];
-0313         else
-0314             model.unconstrained=[model.unconstrained;zeros(numel(metsToAdd),1)];
-0315         end
-0316     else
-0317         if isfield(models{i},'unconstrained')
-0318             model.unconstrained=[zeros(numel(model.mets),1);models{i}.unconstrained(metsToAdd)];
-0319         end
-0320     end
-0321     
-0322     %Only add extra info on new metabolites since it's a little tricky to
-0323     %chose what to keep otherwise. Should change in the future
-0324 
-0325     if ~isempty(metsToAdd)
-0326         if isfield(models{i},'inchis')
-0327             if isfield(model,'inchis')
-0328                 model.inchis=[model.inchis;models{i}.inchis(metsToAdd)];
-0329             else
-0330                 emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
-0331                 emptyInchi(:)={''};
-0332                 model.inchis=[emptyInchi;models{i}.inchis(metsToAdd)];
-0333             end
-0334         else
-0335             if isfield(model,'inchis')
-0336                 emptyInchi=cell(numel(metsToAdd),1);
-0337                 emptyInchi(:)={''};
-0338                 model.inchis=[model.inchis;emptyInchi];
-0339             end
-0340         end
+0097     if ~isempty(conflicting)
+0098         printString=cell(numel(conflicting),1);
+0099         for j=1:numel(conflicting)
+0100             printString{j}=['Old: ' models{i}.rxns{conflicting(j)} ' New: ' models{i}.rxns{conflicting(j)} '_' models{i}.id];
+0101             models{i}.rxns{conflicting(j)}=[models{i}.rxns{conflicting(j)} '_' models{i}.id];
+0102         end
+0103         if supressWarnings==false
+0104             EM=['The following reaction IDs in ' models{i}.id ' are already present in the model and were renamed:'];
+0105             dispEM(EM,false,printString);
+0106             fprintf('\n');
+0107         end
+0108     end
+0109     
+0110     %Add all static stuff
+0111     if any(hasRxnFrom) || (~copyToComps && ~any(hasRxnFrom))
+0112         model.rxnFrom  = [model.rxnFrom;  models{i}.rxnFrom];
+0113     end
+0114     model.rxns     = [model.rxns;     models{i}.rxns];
+0115     model.rxnNames = [model.rxnNames; models{i}.rxnNames];
+0116     model.lb       = [model.lb;       models{i}.lb];
+0117     model.ub       = [model.ub;       models{i}.ub];
+0118     model.c        = [model.c;        models{i}.c];
+0119     model.rev      = [model.rev;      models{i}.rev];
+0120     
+0121     if isfield(models{i},'subSystems')
+0122         if isfield(model,'subSystems')
+0123             model.subSystems=[model.subSystems;models{i}.subSystems];
+0124         else
+0125             emptySubSystem=cell(numel(model.rxns)-numel(models{i}.rxns),1);
+0126             emptySubSystem(:)={''};
+0127             emptySubSystem=cellfun(@(x) cell(0,0),emptySubSystem,'UniformOutput',false);
+0128             model.subSystems=[emptySubSystem;models{i}.subSystems];
+0129         end
+0130     else
+0131         if isfield(model,'subSystems')
+0132             emptySubSystem=cell(numel(models{i}.rxns),1);
+0133             emptySubSystem(:)={''};
+0134             emptySubSystem=cellfun(@(x) cell(0,0),emptySubSystem,'UniformOutput',false);
+0135             model.subSystems=[model.subSystems;emptySubSystem];
+0136         end
+0137     end
+0138     
+0139     if isfield(models{i},'eccodes')
+0140         if isfield(model,'eccodes')
+0141             model.eccodes=[model.eccodes;models{i}.eccodes];
+0142         else
+0143             emptyEC=cell(numel(model.rxns)-numel(models{i}.rxns),1);
+0144             emptyEC(:)={''};
+0145             model.eccodes=[emptyEC;models{i}.eccodes];
+0146         end
+0147     else
+0148         if isfield(model,'eccodes')
+0149             emptyEC=cell(numel(models{i}.rxns),1);
+0150             emptyEC(:)={''};
+0151             model.eccodes=[model.eccodes;emptyEC];
+0152         end
+0153     end
+0154     
+0155     if isfield(models{i},'rxnMiriams')
+0156         if isfield(model,'rxnMiriams')
+0157             model.rxnMiriams=[model.rxnMiriams;models{i}.rxnMiriams];
+0158         else
+0159             model.rxnMiriams=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnMiriams];
+0160         end
+0161     else
+0162         if isfield(model,'rxnMiriams')
+0163             model.rxnMiriams=[model.rxnMiriams;cell(numel(models{i}.rxns),1)];
+0164         end
+0165     end
+0166     
+0167     if isfield(models{i},'rxnNotes')
+0168         if isfield(model,'rxnNotes')
+0169             model.rxnNotes=[model.rxnNotes;models{i}.rxnNotes];
+0170         else
+0171             emptyNotes=cell(numel(model.rxns)-numel(models{i}.rxns),1);
+0172             emptyNotes(:)={''};
+0173             model.rxnNotes=[emptyNotes;models{i}.rxnNotes];
+0174         end
+0175     else
+0176         if isfield(model,'rxnNotes')
+0177             emptyNotes=cell(numel(models{i}.rxns),1);
+0178             emptyNotes(:)={''};
+0179             model.rxnNotes=[model.rxnNotes;emptyNotes];
+0180         end
+0181     end
+0182     
+0183     if isfield(models{i},'rxnReferences')
+0184         if isfield(model,'rxnReferences')
+0185             model.rxnReferences=[model.rxnReferences;models{i}.rxnReferences];
+0186         else
+0187             emptyReferences=cell(numel(model.rxns)-numel(models{i}.rxns),1);
+0188             emptyReferences(:)={''};
+0189             model.rxnReferences=[emptyReferences;models{i}.rxnReferences];
+0190         end
+0191     else
+0192         if isfield(model,'rxnReferences')
+0193             emptyReferences=cell(numel(models{i}.rxns),1);
+0194             emptyReferences(:)={''};
+0195             model.rxnReferences=[model.rxnReferences;emptyReferences];
+0196         end
+0197     end
+0198     
+0199     if isfield(models{i},'rxnConfidenceScores')
+0200         if isfield(model,'rxnConfidenceScores')
+0201             model.rxnConfidenceScores=[model.rxnConfidenceScores;models{i}.rxnConfidenceScores];
+0202         else
+0203             model.rxnConfidenceScores=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnConfidenceScores];
+0204         end
+0205     else
+0206         if isfield(model,'rxnConfidenceScores')
+0207             model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(models{i}.rxns),1)];
+0208         end
+0209     end
+0210 
+0211     if isfield(models{i},'rxnDeltaG')
+0212         if isfield(model,'rxnDeltaG')
+0213             model.rxnDeltaG=[model.rxnDeltaG;models{i}.rxnDeltaG];
+0214         else
+0215             model.rxnDeltaG=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnDeltaG];
+0216         end
+0217     else
+0218         if isfield(model,'rxnDeltaG')
+0219             model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(models{i}.rxns),1)];
+0220         end
+0221     end
+0222     
+0223     if isfield(models{i},'rxnComps')
+0224         if isfield(model,'rxnComps')
+0225             model.rxnComps=[model.rxnComps;models{i}.rxnComps];
+0226         else
+0227             model.rxnComps=[ones(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnComps];
+0228             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');
+0229         end
+0230     else
+0231         if isfield(model,'rxnComps')
+0232             model.rxnComps=[model.rxnComps;ones(numel(models{i}.rxns),1)];
+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     end
+0236     
+0237     if isfield(models{i},'rxnScores')
+0238         if isfield(model,'rxnScores')
+0239             model.rxnScores=[model.rxnScores;models{i}.rxnScores];
+0240         else
+0241             emptyRS=zeros(numel(model.rxns)-numel(models{i}.rxns),1);
+0242             model.rxnScores=[emptyRS;models{i}.rxnScores];
+0243         end
+0244     else
+0245         if isfield(model,'rxnScores')
+0246             emptyRS=zeros(numel(models{i}.rxns),1);
+0247             model.rxnScores=[model.rxnScores;emptyRS];
+0248         end
+0249     end
+0250     
+0251     if isfield(models{i},'pwys')
+0252         if isfield(model,'pwys')
+0253             model.pwys=[model.pwys;models{i}.pwys];
+0254         else
+0255             model.pwys=[cell(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.pwys];
+0256         end
+0257     else
+0258         if isfield(model,'pwys')
+0259             model.pwys=[model.pwys;cell(numel(models{i}.rxns),1)];
+0260         end
+0261     end
+0262 
+0263     if strcmpi(metParam,'metNames')
+0264     %Get the new metabolites from matching the models. Metabolites are said
+0265     %to be the same if they share name and compartment id. This means that
+0266     %metabolite IDs are not taken into account.
+0267         
+0268         oldMetComps=model.comps(model.metComps);
+0269         oldMets=strcat(model.metNames,'[',oldMetComps,']');
+0270         %This is because it makes a '[]' string if no new metabolites
+0271         if ~isempty(models{i}.metNames)
+0272             newMetComps=models{i}.comps(models{i}.metComps);
+0273             newMets=strcat(models{i}.metNames,'[',newMetComps,']');
+0274         else
+0275             newMets={};
+0276             newMetComps={};
+0277         end
+0278         tf=ismember(newMets,oldMets);
+0279         metsToAdd=find(~tf);
+0280 
+0281     end
+0282 
+0283     if strcmpi(metParam,'mets')
+0284     %Get the new metabolites from matching the models. Metabolites are matched by metabolite ID (model.mets).
+0285 
+0286         oldMetComps=model.comps(model.metComps);
+0287         oldMets=model.mets;
+0288     
+0289         if ~isempty(models{i}.mets)
+0290             newMetComps=models{i}.comps(models{i}.metComps);
+0291             newMets=models{i}.mets;
+0292         else
+0293             newMets={};
+0294             newMetComps={};
+0295         end
+0296         tf=ismember(newMets,oldMets);
+0297         metsToAdd=find(~tf);
+0298 
+0299     end
+0300     
+0301     %First add the new metabolites Make sure that there are no conflicting
+0302     %metabolite ids
+0303     conflicting=ismember(models{i}.mets(metsToAdd),model.mets);
+0304     
+0305     conflicting=find(conflicting);
+0306     
+0307     if ~isempty(conflicting)
+0308         printString=cell(numel(conflicting),1);
+0309         for j=1:numel(conflicting)
+0310             printString{j}=['Old: ' models{i}.mets{metsToAdd(conflicting(j))} ' New: ' models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
+0311             models{i}.mets{metsToAdd(conflicting(j))}=[models{i}.mets{metsToAdd(conflicting(j))} '_' models{i}.id];
+0312         end
+0313         if supressWarnings==false
+0314             EM=['The following metabolite IDs in ' models{i}.id ' are already present in the model and were renamed:'];
+0315             dispEM(EM,false,printString);
+0316         end
+0317     end
+0318     
+0319     %Add static info on the metabolites
+0320     if any(hasMetFrom)
+0321         model.metFrom  = [model.metFrom;  models{i}.metFrom(metsToAdd)];
+0322     end
+0323     model.mets     = [model.mets;     models{i}.mets(metsToAdd)];
+0324     model.metNames = [model.metNames; models{i}.metNames(metsToAdd)];
+0325     model.b        = [model.b;        zeros(numel(metsToAdd),size(model.b,2))];
+0326     
+0327     if isfield(model,'unconstrained')
+0328         if isfield(models{i},'unconstrained')
+0329             model.unconstrained=[model.unconstrained;models{i}.unconstrained(metsToAdd)];
+0330         else
+0331             model.unconstrained=[model.unconstrained;zeros(numel(metsToAdd),1)];
+0332         end
+0333     else
+0334         if isfield(models{i},'unconstrained')
+0335             model.unconstrained=[zeros(numel(model.mets),1);models{i}.unconstrained(metsToAdd)];
+0336         end
+0337     end
+0338     
+0339     %Only add extra info on new metabolites since it's a little tricky to
+0340     %chose what to keep otherwise. Should change in the future
 0341 
-0342         if isfield(models{i},'metSmiles')
-0343             if isfield(model,'metSmiles')
-0344                 model.metSmiles=[model.metSmiles;models{i}.metSmiles(metsToAdd)];
-0345             else
-0346                 emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
-0347                 emptyInchi(:)={''};
-0348                 model.metSmiles=[emptyInchi;models{i}.metSmiles(metsToAdd)];
-0349             end
-0350         else
-0351             if isfield(model,'metSmiles')
-0352                 emptyInchi=cell(numel(metsToAdd),1);
-0353                 emptyInchi(:)={''};
-0354                 model.metSmiles=[model.metSmiles;emptyInchi];
-0355             end
-0356         end
-0357         
-0358         if isfield(models{i},'metFormulas')
-0359             if isfield(model,'metFormulas')
-0360                 model.metFormulas=[model.metFormulas;models{i}.metFormulas(metsToAdd)];
-0361             else
-0362                 emptyMetFormulas=cell(numel(model.mets)-numel(metsToAdd),1);
-0363                 emptyMetFormulas(:)={''};
-0364                 model.metFormulas=[emptyMetFormulas;models{i}.metFormulas(metsToAdd)];
-0365             end
-0366         else
-0367             if isfield(model,'metFormulas')
-0368                 emptyMetFormulas=cell(numel(metsToAdd),1);
-0369                 emptyMetFormulas(:)={''};
-0370                 model.metFormulas=[model.metFormulas;emptyMetFormulas];
-0371             end
-0372         end
-0373         
-0374         if isfield(models{i},'metCharges')
-0375             if isfield(model,'metCharges')
-0376                 model.metCharges=[model.metCharges;models{i}.metCharges(metsToAdd)];
-0377             else
-0378                 emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
-0379                 model.metCharges=[emptyMetCharge;models{i}.metCharges(metsToAdd)];
-0380             end
-0381         else
-0382             if isfield(model,'metCharges')
-0383                 emptyMetCharge=nan(numel(metsToAdd),1);
-0384                 model.metCharges=[model.metCharges;emptyMetCharge];
-0385             end
-0386         end
-0387 
-0388         if isfield(models{i},'metDeltaG')
-0389             if isfield(model,'metDeltaG')
-0390                 model.metDeltaG=[model.metDeltaG;models{i}.metDeltaG(metsToAdd)];
-0391             else
-0392                 emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
-0393                 model.metDeltaG=[emptyMetCharge;models{i}.metDeltaG(metsToAdd)];
-0394             end
-0395         else
-0396             if isfield(model,'metDeltaG')
-0397                 emptyMetCharge=nan(numel(metsToAdd),1);
-0398                 model.metDeltaG=[model.metDeltaG;emptyMetCharge];
-0399             end
-0400         end
-0401         
-0402         if isfield(models{i},'metMiriams')
-0403             if isfield(model,'metMiriams')
-0404                 model.metMiriams=[model.metMiriams;models{i}.metMiriams(metsToAdd)];
-0405             else
-0406                 emptyMetMiriam=cell(numel(model.mets)-numel(metsToAdd),1);
-0407                 model.metMiriams=[emptyMetMiriam;models{i}.metMiriams(metsToAdd)];
-0408             end
-0409         else
-0410             if isfield(model,'metMiriams')
-0411                 emptyMetMiriam=cell(numel(metsToAdd),1);
-0412                 model.metMiriams=[model.metMiriams;emptyMetMiriam];
-0413             end
-0414         end
-0415     end
-0416     
-0417     %Add if there are any new compartments and add those. This can change
-0418     %the order of compartments and the corresponding indexes in
-0419     %model.metComps.
-0420     
-0421     %Find overlapping and new compartments
-0422     [overlap, oldIDs]=ismember(models{i}.comps,model.comps);
-0423     overlap=find(overlap);
-0424     
-0425     %Add the new compartments if any
-0426     if numel(overlap)~=numel(models{i}.compNames)
-0427         compIndexes=oldIDs==0;
-0428         
-0429         %Make sure that there are no conflicting compartment ids
-0430         [~, conflicting]=ismember(models{i}.compNames(compIndexes),model.compNames);
-0431         if any(conflicting)
-0432             EM=['The following compartment IDs in ' models{i}.id ' are already present in the model but with another name. They have to be renamed'];
-0433             dispEM(EM,true,model.comps(conflicting));
-0434         end
-0435         
-0436         %It's ok to add duplicate name, but not duplicate IDs
-0437         model.compNames=[model.compNames; models{i}.compNames(compIndexes)];
-0438         model.comps=[model.comps; models{i}.comps(compIndexes)];
-0439         if isfield(model,'compOutside')
-0440             if isfield(models{i},'compOutside')
-0441                 model.compOutside=[model.compOutside; models{i}.compOutside(compIndexes)];
-0442             else
-0443                 %This is if not all models have the field
-0444                 padding=cell(sum(compIndexes),1);
-0445                 padding(:)={''};
-0446                 model.compOutside=[model.compOutside;padding];
-0447             end
-0448         end
-0449         if isfield(model,'compMiriams')
-0450             if isfield(models{i},'compMiriams')
-0451                 model.compMiriams=[model.compMiriams; models{i}.compMiriams(compIndexes)];
-0452             else
-0453                 %This is if not all models have the field
-0454                 model.compMiriams=[model.compMiriams;cell(sum(compIndexes),1)];
-0455             end
-0456         end
-0457     end
-0458     
-0459     %Only add new comp info on the un-matched metabolites since the old
-0460     %ones will be mapped to the existing list anyways
-0461     [I, J]=ismember(newMetComps(metsToAdd),model.comps);
-0462     %Just a check
-0463     if ~all(I)
-0464         EM='There was an unexpected error in matching compartments';
-0465         dispEM(EM);
-0466     end
-0467     model.metComps=[model.metComps;J];
-0468      
-0469     %Create the new stoichiometric matrix
-0470     model.S=[model.S;sparse(numel(metsToAdd),size(model.S,2))];
-0471     
-0472 
-0473     if strcmpi(metParam,'metNames')
-0474         %Rematch metabolite names. Not the most clever way to do it maybe
-0475         allMets=strcat(model.metNames,'[',model.comps(model.metComps),']');
-0476         [~, J]=ismember(newMets,allMets);
-0477     end
-0478 
-0479     if strcmpi(metParam,'mets')
-0480         %Rematch metabolite by IDs and add unique new metabolites
-0481         allMets=model.mets;
-0482         uniqueNewMets = setdiff(newMets,oldMets);
-0483         allMets(end+1:end+numel(uniqueNewMets)) = uniqueNewMets;
-0484         [~, J]=ismember(newMets,allMets);
-0485     end
-0486 
-0487     %Update the stoichiometric matrix for the model to add
-0488     newS=sparse(numel(model.mets),numel(models{i}.rxns));
-0489     newS(J,:)=models{i}.S;
-0490     model.S=[model.S newS];
-0491 
-0492     
-0493     %Now add new genes
-0494     if isfield(models{i},'genes')
-0495         if ~isfield(model,'genes')
-0496             %If there was no gene info before
-0497             model.genes      = models{i}.genes;
-0498             model.rxnGeneMat = [sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
-0499             emptyGene        = repmat({''},numel(model.rxns),1);
-0500             model.grRules    = [emptyGene;models{i}.grRules];
-0501             model.geneFrom   = models{i}.geneFrom;
-0502             
-0503             if isfield(models{i},'geneShortNames')
-0504                 model.geneShortNames=models{i}.geneShortNames;
-0505             end
-0506 
-0507             if isfield(models{i},'proteins')
-0508                 model.proteins=models{i}.proteins;
-0509             end
-0510 
-0511             if isfield(models{i},'geneMiriams')
-0512                 model.geneMiriams=models{i}.geneMiriams;
-0513             end
-0514             
-0515             if isfield(models{i},'geneComps')
-0516                 model.geneComps=models{i}.geneComps;
-0517             end
-0518         else
-0519             %If gene info should be merged
-0520             a=ismember(models{i}.genes,model.genes);
+0342     if ~isempty(metsToAdd)
+0343         if isfield(models{i},'inchis')
+0344             if isfield(model,'inchis')
+0345                 model.inchis=[model.inchis;models{i}.inchis(metsToAdd)];
+0346             else
+0347                 emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
+0348                 emptyInchi(:)={''};
+0349                 model.inchis=[emptyInchi;models{i}.inchis(metsToAdd)];
+0350             end
+0351         else
+0352             if isfield(model,'inchis')
+0353                 emptyInchi=cell(numel(metsToAdd),1);
+0354                 emptyInchi(:)={''};
+0355                 model.inchis=[model.inchis;emptyInchi];
+0356             end
+0357         end
+0358 
+0359         if isfield(models{i},'metSmiles')
+0360             if isfield(model,'metSmiles')
+0361                 model.metSmiles=[model.metSmiles;models{i}.metSmiles(metsToAdd)];
+0362             else
+0363                 emptyInchi=cell(numel(model.mets)-numel(metsToAdd),1);
+0364                 emptyInchi(:)={''};
+0365                 model.metSmiles=[emptyInchi;models{i}.metSmiles(metsToAdd)];
+0366             end
+0367         else
+0368             if isfield(model,'metSmiles')
+0369                 emptyInchi=cell(numel(metsToAdd),1);
+0370                 emptyInchi(:)={''};
+0371                 model.metSmiles=[model.metSmiles;emptyInchi];
+0372             end
+0373         end
+0374         
+0375         if isfield(models{i},'metFormulas')
+0376             if isfield(model,'metFormulas')
+0377                 model.metFormulas=[model.metFormulas;models{i}.metFormulas(metsToAdd)];
+0378             else
+0379                 emptyMetFormulas=cell(numel(model.mets)-numel(metsToAdd),1);
+0380                 emptyMetFormulas(:)={''};
+0381                 model.metFormulas=[emptyMetFormulas;models{i}.metFormulas(metsToAdd)];
+0382             end
+0383         else
+0384             if isfield(model,'metFormulas')
+0385                 emptyMetFormulas=cell(numel(metsToAdd),1);
+0386                 emptyMetFormulas(:)={''};
+0387                 model.metFormulas=[model.metFormulas;emptyMetFormulas];
+0388             end
+0389         end
+0390         
+0391         if isfield(models{i},'metCharges')
+0392             if isfield(model,'metCharges')
+0393                 model.metCharges=[model.metCharges;models{i}.metCharges(metsToAdd)];
+0394             else
+0395                 emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
+0396                 model.metCharges=[emptyMetCharge;models{i}.metCharges(metsToAdd)];
+0397             end
+0398         else
+0399             if isfield(model,'metCharges')
+0400                 emptyMetCharge=nan(numel(metsToAdd),1);
+0401                 model.metCharges=[model.metCharges;emptyMetCharge];
+0402             end
+0403         end
+0404 
+0405         if isfield(models{i},'metDeltaG')
+0406             if isfield(model,'metDeltaG')
+0407                 model.metDeltaG=[model.metDeltaG;models{i}.metDeltaG(metsToAdd)];
+0408             else
+0409                 emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
+0410                 model.metDeltaG=[emptyMetCharge;models{i}.metDeltaG(metsToAdd)];
+0411             end
+0412         else
+0413             if isfield(model,'metDeltaG')
+0414                 emptyMetCharge=nan(numel(metsToAdd),1);
+0415                 model.metDeltaG=[model.metDeltaG;emptyMetCharge];
+0416             end
+0417         end
+0418         
+0419         if isfield(models{i},'metMiriams')
+0420             if isfield(model,'metMiriams')
+0421                 model.metMiriams=[model.metMiriams;models{i}.metMiriams(metsToAdd)];
+0422             else
+0423                 emptyMetMiriam=cell(numel(model.mets)-numel(metsToAdd),1);
+0424                 model.metMiriams=[emptyMetMiriam;models{i}.metMiriams(metsToAdd)];
+0425             end
+0426         else
+0427             if isfield(model,'metMiriams')
+0428                 emptyMetMiriam=cell(numel(metsToAdd),1);
+0429                 model.metMiriams=[model.metMiriams;emptyMetMiriam];
+0430             end
+0431         end
+0432     end
+0433     
+0434     %Add if there are any new compartments and add those. This can change
+0435     %the order of compartments and the corresponding indexes in
+0436     %model.metComps.
+0437     
+0438     %Find overlapping and new compartments
+0439     [overlap, oldIDs]=ismember(models{i}.comps,model.comps);
+0440     overlap=find(overlap);
+0441     
+0442     %Add the new compartments if any
+0443     if numel(overlap)~=numel(models{i}.compNames)
+0444         compIndexes=oldIDs==0;
+0445         
+0446         %Make sure that there are no conflicting compartment ids
+0447         [~, conflicting]=ismember(models{i}.compNames(compIndexes),model.compNames);
+0448         if any(conflicting)
+0449             EM=['The following compartment IDs in ' models{i}.id ' are already present in the model but with another name. They have to be renamed'];
+0450             dispEM(EM,true,model.comps(conflicting));
+0451         end
+0452         
+0453         %It's ok to add duplicate name, but not duplicate IDs
+0454         model.compNames=[model.compNames; models{i}.compNames(compIndexes)];
+0455         model.comps=[model.comps; models{i}.comps(compIndexes)];
+0456         if isfield(model,'compOutside')
+0457             if isfield(models{i},'compOutside')
+0458                 model.compOutside=[model.compOutside; models{i}.compOutside(compIndexes)];
+0459             else
+0460                 %This is if not all models have the field
+0461                 padding=cell(sum(compIndexes),1);
+0462                 padding(:)={''};
+0463                 model.compOutside=[model.compOutside;padding];
+0464             end
+0465         end
+0466         if isfield(model,'compMiriams')
+0467             if isfield(models{i},'compMiriams')
+0468                 model.compMiriams=[model.compMiriams; models{i}.compMiriams(compIndexes)];
+0469             else
+0470                 %This is if not all models have the field
+0471                 model.compMiriams=[model.compMiriams;cell(sum(compIndexes),1)];
+0472             end
+0473         end
+0474     end
+0475     
+0476     %Only add new comp info on the un-matched metabolites since the old
+0477     %ones will be mapped to the existing list anyways
+0478     [I, J]=ismember(newMetComps(metsToAdd),model.comps);
+0479     %Just a check
+0480     if ~all(I)
+0481         EM='There was an unexpected error in matching compartments';
+0482         dispEM(EM);
+0483     end
+0484     model.metComps=[model.metComps;J];
+0485      
+0486     %Create the new stoichiometric matrix
+0487     model.S=[model.S;sparse(numel(metsToAdd),size(model.S,2))];
+0488     
+0489 
+0490     if strcmpi(metParam,'metNames')
+0491         %Rematch metabolite names. Not the most clever way to do it maybe
+0492         allMets=strcat(model.metNames,'[',model.comps(model.metComps),']');
+0493         [~, J]=ismember(newMets,allMets);
+0494     end
+0495 
+0496     if strcmpi(metParam,'mets')
+0497         %Rematch metabolite by IDs and add unique new metabolites
+0498         allMets=model.mets;
+0499         uniqueNewMets = setdiff(newMets,oldMets);
+0500         allMets(end+1:end+numel(uniqueNewMets)) = uniqueNewMets;
+0501         [~, J]=ismember(newMets,allMets);
+0502     end
+0503 
+0504     %Update the stoichiometric matrix for the model to add
+0505     newS=sparse(numel(model.mets),numel(models{i}.rxns));
+0506     newS(J,:)=models{i}.S;
+0507     model.S=[model.S newS];
+0508 
+0509     
+0510     %Now add new genes
+0511     if isfield(models{i},'genes')
+0512         if ~isfield(model,'genes')
+0513             %If there was no gene info before
+0514             model.genes      = models{i}.genes;
+0515             model.rxnGeneMat = [sparse(numel(model.rxns),numel(models{i}.genes));models{i}.rxnGeneMat];
+0516             emptyGene        = repmat({''},numel(model.rxns),1);
+0517             model.grRules    = [emptyGene;models{i}.grRules];
+0518             if any(hasGeneFrom)
+0519                 model.geneFrom   = models{i}.geneFrom;
+0520             end
 0521             
-0522             genesToAdd=find(~a);
-0523             
-0524             %Only add extra gene info on new genes. This might not be
-0525             %correct and should be changed later...
-0526             if ~isempty(genesToAdd)
-0527                 model.genes      = [model.genes;        models{i}.genes(genesToAdd)];
-0528                 model.geneFrom   = [model.geneFrom;     models{i}.geneFrom(genesToAdd)];
-0529                 model.rxnGeneMat = [model.rxnGeneMat    sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
-0530                 
-0531                 if isfield(models{i},'geneShortNames')
-0532                     if isfield(model,'geneShortNames')
-0533                         model.geneShortNames=[model.geneShortNames;models{i}.geneShortNames(genesToAdd)];
-0534                     else
-0535                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
-0536                         emptyGeneSN(:)={''};
-0537                         model.geneShortNames=[emptyGeneSN;models{i}.geneShortNames(genesToAdd)];
-0538                     end
-0539                 else
-0540                     if isfield(model,'geneShortNames')
-0541                         emptyGeneSN=cell(numel(genesToAdd),1);
-0542                         emptyGeneSN(:)={''};
-0543                         model.geneShortNames=[model.geneShortNames;emptyGeneSN];
-0544                     end
-0545                 end
-0546 
-0547                 if isfield(models{i},'proteins')
-0548                     if isfield(model,'proteins')
-0549                         model.proteins=[model.proteins;models{i}.proteins(genesToAdd)];
-0550                     else
-0551                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
-0552                         emptyGeneSN(:)={''};
-0553                         model.proteins=[emptyGeneSN;models{i}.proteins(genesToAdd)];
-0554                     end
-0555                 else
-0556                     if isfield(model,'proteins')
-0557                         emptyGeneSN=cell(numel(genesToAdd),1);
-0558                         emptyGeneSN(:)={''};
-0559                         model.proteins=[model.proteins;emptyGeneSN];
-0560                     end
-0561                 end
-0562 
-0563                 if isfield(models{i},'geneMiriams')
-0564                     if isfield(model,'geneMiriams')
-0565                         model.geneMiriams=[model.geneMiriams;models{i}.geneMiriams(genesToAdd)];
-0566                     else
-0567                         emptyGeneMir=cell(numel(model.genes)-numel(genesToAdd),1);
-0568                         model.geneMiriams=[emptyGeneMir;models{i}.geneMiriams(genesToAdd)];
-0569                     end
-0570                 else
-0571                     if isfield(model,'geneMiriams')
-0572                         emptyGeneMir=cell(numel(genesToAdd),1);
-0573                         model.geneMiriams=[model.geneMiriams;emptyGeneMir];
-0574                     end
-0575                 end
-0576                 
-0577                 if isfield(models{i},'geneComps')
-0578                     if isfield(model,'geneComps')
-0579                         model.geneComps=[model.geneComps;models{i}.geneComps(genesToAdd)];
-0580                     else
-0581                         emptyGeneMir=ones(numel(model.genes)-numel(genesToAdd),1);
-0582                         model.geneComps=[emptyGeneMir;models{i}.geneComps(genesToAdd)];
-0583                         EM='Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment';
-0584                         dispEM(EM,false);
-0585                     end
-0586                 else
-0587                     if isfield(model,'geneComps')
-0588                         emptyGeneMir=ones(numel(genesToAdd),1);
-0589                         model.geneComps=[model.geneComps;emptyGeneMir];
-0590                         EM='Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment';
-0591                         dispEM(EM,false);
-0592                     end
-0593                 end
-0594             end
-0595             
-0596             %Remap the genes from the new model. The same thing as with
-0597             %mets; this is a wasteful way to do it but I don't care right
-0598             %now
-0599             a = ismember(models{i}.genes,model.genes);
-0600             
-0601             %Just a check
-0602             if ~all(a)
-0603                 EM='There was an unexpected error in matching genes';
-0604                 dispEM(EM);
-0605             end
-0606             model.grRules=[model.grRules;models{i}.grRules];
-0607         end
-0608     else
-0609         %Add empty gene associations
-0610         if isfield(model,'genes')
-0611             emptyGene=cell(numel(models{i}.rxns),1);
-0612             emptyGene(:)={''};
-0613             model.grRules=[model.grRules;emptyGene];
-0614         end
-0615     end
-0616 end
-0617 %Fix grRules and reconstruct rxnGeneMat
-0618 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
-0619 model.grRules = grRules;
-0620 model.rxnGeneMat = rxnGeneMat;
-0621 end
    
+0522             if isfield(models{i},'geneShortNames')
+0523                 model.geneShortNames=models{i}.geneShortNames;
+0524             end
+0525 
+0526             if isfield(models{i},'proteins')
+0527                 model.proteins=models{i}.proteins;
+0528             end
+0529 
+0530             if isfield(models{i},'geneMiriams')
+0531                 model.geneMiriams=models{i}.geneMiriams;
+0532             end
+0533             
+0534             if isfield(models{i},'geneComps')
+0535                 model.geneComps=models{i}.geneComps;
+0536             end
+0537         else
+0538             %If gene info should be merged
+0539             a=ismember(models{i}.genes,model.genes);
+0540             
+0541             genesToAdd=find(~a);
+0542             
+0543             %Only add extra gene info on new genes. This might not be
+0544             %correct and should be changed later...
+0545             if ~isempty(genesToAdd)
+0546                 model.genes      = [model.genes;        models{i}.genes(genesToAdd)];
+0547                 model.geneFrom   = [model.geneFrom;     models{i}.geneFrom(genesToAdd)];
+0548                 model.rxnGeneMat = [model.rxnGeneMat    sparse(size(model.rxnGeneMat,1),numel(genesToAdd))];
+0549                 
+0550                 if isfield(models{i},'geneShortNames')
+0551                     if isfield(model,'geneShortNames')
+0552                         model.geneShortNames=[model.geneShortNames;models{i}.geneShortNames(genesToAdd)];
+0553                     else
+0554                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
+0555                         emptyGeneSN(:)={''};
+0556                         model.geneShortNames=[emptyGeneSN;models{i}.geneShortNames(genesToAdd)];
+0557                     end
+0558                 else
+0559                     if isfield(model,'geneShortNames')
+0560                         emptyGeneSN=cell(numel(genesToAdd),1);
+0561                         emptyGeneSN(:)={''};
+0562                         model.geneShortNames=[model.geneShortNames;emptyGeneSN];
+0563                     end
+0564                 end
+0565 
+0566                 if isfield(models{i},'proteins')
+0567                     if isfield(model,'proteins')
+0568                         model.proteins=[model.proteins;models{i}.proteins(genesToAdd)];
+0569                     else
+0570                         emptyGeneSN=cell(numel(model.genes)-numel(genesToAdd),1);
+0571                         emptyGeneSN(:)={''};
+0572                         model.proteins=[emptyGeneSN;models{i}.proteins(genesToAdd)];
+0573                     end
+0574                 else
+0575                     if isfield(model,'proteins')
+0576                         emptyGeneSN=cell(numel(genesToAdd),1);
+0577                         emptyGeneSN(:)={''};
+0578                         model.proteins=[model.proteins;emptyGeneSN];
+0579                     end
+0580                 end
+0581 
+0582                 if isfield(models{i},'geneMiriams')
+0583                     if isfield(model,'geneMiriams')
+0584                         model.geneMiriams=[model.geneMiriams;models{i}.geneMiriams(genesToAdd)];
+0585                     else
+0586                         emptyGeneMir=cell(numel(model.genes)-numel(genesToAdd),1);
+0587                         model.geneMiriams=[emptyGeneMir;models{i}.geneMiriams(genesToAdd)];
+0588                     end
+0589                 else
+0590                     if isfield(model,'geneMiriams')
+0591                         emptyGeneMir=cell(numel(genesToAdd),1);
+0592                         model.geneMiriams=[model.geneMiriams;emptyGeneMir];
+0593                     end
+0594                 end
+0595                 
+0596                 if isfield(models{i},'geneComps')
+0597                     if isfield(model,'geneComps')
+0598                         model.geneComps=[model.geneComps;models{i}.geneComps(genesToAdd)];
+0599                     else
+0600                         emptyGeneMir=ones(numel(model.genes)-numel(genesToAdd),1);
+0601                         model.geneComps=[emptyGeneMir;models{i}.geneComps(genesToAdd)];
+0602                         EM='Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment';
+0603                         dispEM(EM,false);
+0604                     end
+0605                 else
+0606                     if isfield(model,'geneComps')
+0607                         emptyGeneMir=ones(numel(genesToAdd),1);
+0608                         model.geneComps=[model.geneComps;emptyGeneMir];
+0609                         EM='Adding genes with compartment information to a model without such information. All existing genes will be assigned to the first compartment';
+0610                         dispEM(EM,false);
+0611                     end
+0612                 end
+0613             end
+0614             
+0615             %Remap the genes from the new model. The same thing as with
+0616             %mets; this is a wasteful way to do it but I don't care right
+0617             %now
+0618             a = ismember(models{i}.genes,model.genes);
+0619             
+0620             %Just a check
+0621             if ~all(a)
+0622                 EM='There was an unexpected error in matching genes';
+0623                 dispEM(EM);
+0624             end
+0625             model.grRules=[model.grRules;models{i}.grRules];
+0626         end
+0627     else
+0628         %Add empty gene associations
+0629         if isfield(model,'genes')
+0630             emptyGene=cell(numel(models{i}.rxns),1);
+0631             emptyGene(:)={''};
+0632             model.grRules=[model.grRules;emptyGene];
+0633         end
+0634     end
+0635 end
+0636 %Fix grRules and reconstruct rxnGeneMat
+0637 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
+0638 model.grRules = grRules;
+0639 model.rxnGeneMat = rxnGeneMat;
+0640 end
    
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file
diff --git a/doc/core/standardizeGrRules.html b/doc/core/standardizeGrRules.html
index b9cb023f..b58d66bd 100644
--- a/doc/core/standardizeGrRules.html
+++ b/doc/core/standardizeGrRules.html
@@ -56,7 +56,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 </ul>
 This function is called by:
 <ul style=
-
  • addRxns	addRxns
  • changeGrRules	changeGrRules
  • contractModel	contractModel
  • expandModel	expandModel
  • getModelFromHomology	getModelFromHomology
  • mergeModels	mergeModels
  • predictLocalization	predictLocalization
  • removeGenes	removeGenes
    • 
+
  • addRxns	addRxns
  • changeGrRules	changeGrRules
  • contractModel	contractModel
  • expandModel	expandModel
  • getModelFromHomology	getModelFromHomology
  • mergeModels	mergeModels
  • predictLocalization	predictLocalization
  • removeGenes	removeGenes
    • 
 
 
 
    SUBFUNCTIONS ^
    
diff --git a/doc/index.html b/doc/index.html
index 25c6ef92..faf716ce 100644
--- a/doc/index.html
+++ b/doc/index.html
@@ -3,8 +3,8 @@
 
 
   Matlab Index
-  
-  
+  
+  
   
   
   
@@ -15,63 +15,64 @@
 
    Matlab Index
    
 
    Matlab Directories
    
 
+
  • INIT
  • core
  • external
  • external\kegg
  • external\metacyc
  • hpa
  • installation
  • io
  • legacy\core
  • legacy\external
  • pathway
  • plotting
  • solver
  • struct_conversion
  • testing\manual_tests
  • testing\unit_tests
  • tutorial
  • utils
    • 
 
    Matlab Files found in these Directories
    
 
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
-				
    
 		
    FSEOF		editMiriam		getModelFromMetaCyc		qMOMA	
    	
    INITStepDesc		expandModel		getObjectiveString		randomSampling	
    	
    ManualINITTests		exportForGit		getPathwayDimensions		ravenCobraWrapper	
    	
    SBMLFromExcel		exportModel		getPhylDist		readYAMLmodel	
    	
    addExchangeRxns		exportModelToSIF		getRxnsFromKEGG		removeBadRxns	
    	
    addGenesRaven		exportToExcelFormat		getRxnsFromMetaCyc		removeGenes	
    	
    addIdentifierPrefix		exportToTabDelimited		getRxnsInComp		removeIdentifierPrefix	
    	
    addJavaPaths		extractMiriam		getToolboxVersion		removeLowScoreGenes	
    	
    addMets		fillGaps		getTransportRxns		removeMets	
    	
    addRavenToUserPath		fillGapsLargeTests		getWSLpath		removeRavenFromPath	
    	
    addRxns		fillGapsSmallTests		getWoLFScores		removeReactions	
    	
    addRxnsGenesMets		findGeneDeletions		groupRxnScores		replaceMets	
    	
    addSpontaneousRxns		findRAVENroot		guessComposition		reporterMetabolites	
    	
    addTransport		fitParameters		haveFlux		rescaleModelForINIT	
    	
    analyzeSampling		fitTasks		hmmerTests		reverseRxns	
    	
    blastPlusTests		followChanged		importExcelModel		runDynamicFBA	
    	
    buildEquation		followFluxes		importExportTests		runINIT	
    	
    canConsume		ftINIT		importModel		runPhenotypePhasePlane	
    	
    canProduce		ftINITFillGaps		linkMetaCycKEGGRxns		runProductionEnvelope	
    	
    cdhitTests		ftINITFillGapsForAllTasks		loadSheet		runRobustnessAnalysis	
    	
    changeGeneAssoc		ftINITFillGapsMILP		loadWorkbook		runSimpleOptKnock	
    	
    changeGrRules		ftINITInternalAlg		mafftTests		scoreComplexModel	
    	
    changeRxns		gapReport		makeFakeBlastStructure		scoreModel	
    	
    checkFileExistence		generateNewIds		makeSomething		setColorToMapRxns	
    	
    checkFunctionUniqueness		getAllRxnsFromGenes		mapCompartments		setExchangeBounds	
    	
    checkInstallation		getAllSubGraphs		mapPathwayRxnNames		setOmicDataToRxns	
    	
    checkModelStruct		getAllowedBounds		markPathwayWithExpression		setParam	
    	
    checkProduction		getBlast		markPathwayWithFluxes		setRavenSolver	
    	
    checkRxn		getBlastFromExcel		mergeCompartments		setTitle	
    	
    checkSolution		getColorCodes		mergeLinear		simplifyModel	
    	
    checkTasks		getDiamond		mergeModels		solveLP	
    	
    checkTasksTests		getElementalBalance		miriamTests		solveQP	
    	
    cleanSheet		getEnzymesFromMetaCyc		modelAbilitiesTests		solverTests	
    	
    closeModel		getEssentialRxns		modelConversionTests		sortIdentifiers	
    	
    colorPathway		getExchangeRxns		modelCurationTests		sortModel	
    	
    colorSubsystem		getExprForRxnScore		modelSortingTests		standardizeGrRules	
    	
    combineMetaCycKEGGModels		getExpressionStructure		optimizeProb		standardizeModelFieldOrder	
    	
    compareMultipleModels		getFluxZ		parallelPoolRAVEN		startup	
    	
    compareRxnsGenesMetsComps		getFullPath		parseFormulas		tinitTests	
    	
    constructEquations		getGenesFromGrRules		parseHPA		trimPathway	
    	
    constructMultiFasta		getGenesFromKEGG		parseHPArna		tutorial1	
    	
    constructPathwayFromCelldesigner		getINITModel		parseRxnEqu		tutorial2	
    	
    constructS		getINITSteps		parseScores		tutorial2_solutions	
    	
    consumeSomething		getIndexes		parseTaskList		tutorial3	
    	
    contractModel		getKEGGModelForOrganism		permuteModel		tutorial3_solutions	
    	
    convertCharArray		getMD5Hash		plotAdditionalInfo		tutorial4	
    	
    convertToIrrev		getMetaCycModelForOrganism		plotLabels		tutorial4_solutions	
    	
    copyToComps		getMetsFromKEGG		predictLocalization		tutorial5	
    	
    deleteUnusedGenes		getMetsFromMetaCyc		prepINITModel		tutorial6	
    	
    diamondTests		getMetsInComp		printFluxes		updateDocumentation	
    	
    dispEM		getMinNrFluxes		printModel		writeSheet	
    	
    drawMap		getModelFromHomology		printModelStats		writeYAMLmodel	
    	
    drawPathway		getModelFromKEGG		printOrange			
    
+				FSEOF		emptyOrLogicalScalar		getModelFromKEGG		qMOMA		
+				INITStepDesc		emptyOrTextOrCellOfText		getModelFromMetaCyc		randomSampling		
+				ManualINITTests		emptyOrTextScalar		getObjectiveString		ravenCobraWrapper		
+				SBMLFromExcel		expandModel		getPathwayDimensions		readYAMLmodel		
+				addExchangeRxns		exportForGit		getPhylDist		removeBadRxns		
+				addGenesRaven		exportModel		getRxnsFromKEGG		removeGenes		
+				addIdentifierPrefix		exportModelToSIF		getRxnsFromMetaCyc		removeIdentifierPrefix		
+				addJavaPaths		exportToExcelFormat		getRxnsInComp		removeLowScoreGenes		
+				addMets		exportToTabDelimited		getToolboxVersion		removeMets		
+				addRavenToUserPath		extractMiriam		getTransportRxns		removeRavenFromPath		
+				addRxns		fillGaps		getWSLpath		removeReactions		
+				addRxnsGenesMets		fillGapsLargeTests		getWoLFScores		replaceMets		
+				addSpontaneousRxns		fillGapsSmallTests		groupRxnScores		reporterMetabolites		
+				addTransport		findGeneDeletions		guessComposition		rescaleModelForINIT		
+				analyzeSampling		findRAVENroot		haveFlux		reverseRxns		
+				blastPlusTests		fitParameters		hmmerTests		runDynamicFBA		
+				buildEquation		fitTasks		importExcelModel		runINIT		
+				canConsume		followChanged		importExportTests		runPhenotypePhasePlane		
+				canProduce		followFluxes		importModel		runProductionEnvelope		
+				cdhitTests		ftINIT		linkMetaCycKEGGRxns		runRobustnessAnalysis		
+				changeGeneAssoc		ftINITFillGaps		loadSheet		runSimpleOptKnock		
+				changeGrRules		ftINITFillGapsForAllTasks		loadWorkbook		scoreComplexModel		
+				changeRxns		ftINITFillGapsMILP		mafftTests		scoreModel		
+				checkFileExistence		ftINITInternalAlg		makeFakeBlastStructure		setColorToMapRxns		
+				checkFunctionUniqueness		gapReport		makeSomething		setExchangeBounds		
+				checkInstallation		generateNewIds		mapCompartments		setOmicDataToRxns		
+				checkModelStruct		getAllRxnsFromGenes		mapPathwayRxnNames		setParam		
+				checkProduction		getAllSubGraphs		markPathwayWithExpression		setRavenSolver		
+				checkRxn		getAllowedBounds		markPathwayWithFluxes		setTitle		
+				checkSolution		getBlast		mergeCompartments		simplifyModel		
+				checkTasks		getBlastFromExcel		mergeLinear		solveLP		
+				checkTasksTests		getColorCodes		mergeModels		solveQP		
+				cleanSheet		getDiamond		miriamTests		solverTests		
+				closeModel		getElementalBalance		modelAbilitiesTests		sortIdentifiers		
+				colorPathway		getEnzymesFromMetaCyc		modelConversionTests		sortModel		
+				colorSubsystem		getEssentialRxns		modelCurationTests		standardizeGrRules		
+				combineMetaCycKEGGModels		getExchangeRxns		modelSortingTests		standardizeModelFieldOrder		
+				compareMultipleModels		getExprForRxnScore		optimizeProb		startup		
+				compareRxnsGenesMetsComps		getExpressionStructure		parallelPoolRAVEN		tinitTests		
+				constructEquations		getFluxZ		parseFormulas		trimPathway		
+				constructMultiFasta		getFullPath		parseHPA		tutorial1		
+				constructPathwayFromCelldesigner		getGenesFromGrRules		parseHPArna		tutorial2		
+				constructS		getGenesFromKEGG		parseRxnEqu		tutorial2_solutions		
+				consumeSomething		getINITModel		parseScores		tutorial3		
+				contractModel		getINITSteps		parseTaskList		tutorial3_solutions		
+				convertCharArray		getIndexes		permuteModel		tutorial4		
+				convertToIrrev		getKEGGModelForOrganism		plotAdditionalInfo		tutorial4_solutions		
+				copyToComps		getMD5Hash		plotLabels		tutorial5		
+				deleteUnusedGenes		getMetaCycModelForOrganism		predictLocalization		tutorial6		
+				diamondTests		getMetsFromKEGG		prepINITModel		updateDocumentation		
+				dispEM		getMetsFromMetaCyc		printFluxes		writeSheet		
+				drawMap		getMetsInComp		printModel		writeYAMLmodel		
+				drawPathway		getMinNrFluxes		printModelStats				
+				editMiriam		getModelFromHomology		printOrange			
 
 
 
    Generated by m2html © 2005
    
diff --git a/doc/testing/unit_tests/importExportTests.html b/doc/testing/unit_tests/importExportTests.html
index 6e32b08f..60dc50cf 100644
--- a/doc/testing/unit_tests/importExportTests.html
+++ b/doc/testing/unit_tests/importExportTests.html
@@ -83,7 +83,7 @@ 
    SOURCE CODE ^function testExcelExport(testCase)
 0037 sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
 0038 load(fullfile(sourceDir,'testing','unit_tests','test_data','ecoli_textbook.mat'), 'model');
-0039 exportToExcelFormat(model,fullfile(sourceDir,'testing','unit_tests','test_data','_test.xlsx'));
+0039 evalc('exportToExcelFormat(model,fullfile(sourceDir,''testing'',''unit_tests'',''test_data'',''_test.xlsx''))');
 0040 %File will not be exactly equal as it contains the current date and time,
 0041 %so md5 or similar would not work. Just check whether file is reasonably
 0042 %sized.
diff --git a/doc/tutorial/index.html b/doc/tutorial/index.html
index b1ccd2cc..efef9dc9 100644
--- a/doc/tutorial/index.html
+++ b/doc/tutorial/index.html
@@ -24,7 +24,9 @@ 
    Matlab files in this directory:
    
 
    Other Matlab-specific files in this directory:
    
 
      
 
    • empty.mat
    • iMK1208+suppInfo.mat
    • pathway.mat
    • pcPathway.mat
    
-
+
    Subsequent directories:
    
+
      
+
    • struct_conversion
    
 
 
    Generated by m2html © 2005
    
 
diff --git a/doc/utils/emptyOrLogicalScalar.html b/doc/utils/emptyOrLogicalScalar.html
new file mode 100644
index 00000000..a7371337
--- /dev/null
+++ b/doc/utils/emptyOrLogicalScalar.html
@@ -0,0 +1,53 @@
+
+
+
+  Description of emptyOrLogicalScalar
+  
+  
+  
+  
+  
+  
+
+
+
+
    Home >  utils > emptyOrLogicalScalar.m
    
+
+
+
+
    emptyOrLogicalScalar
+
    
+
+
    PURPOSE ^
    
+
    
+
+
    SYNOPSIS ^
    
+
    function emptyOrLogical(x) 
    
+
+
    DESCRIPTION ^
    
+
    
+
+
+
    CROSS-REFERENCE INFORMATION ^
    
+This function calls:
+
      
+
    
+This function is called by:
+
      
+
    
+
+
+
+
+
    SOURCE CODE ^
    
+
    0001 function emptyOrLogical(x)
+0002     if isempty(x)
+0003         return;
+0004     end
+0005     validateattributes(x, {'logical'}, {'scalar'});
+0006 end
    
+
    Generated by m2html © 2005
    
+
+
\ No newline at end of file
diff --git a/doc/utils/emptyOrTextOrCellOfText.html b/doc/utils/emptyOrTextOrCellOfText.html
new file mode 100644
index 00000000..72298c9a
--- /dev/null
+++ b/doc/utils/emptyOrTextOrCellOfText.html
@@ -0,0 +1,90 @@
+
+
+
+  Description of emptyOrTextOrCellOfText
+  
+  
+  
+  
+  
+  
+
+
+
+
    Home >  utils > emptyOrTextOrCellOfText.m
    
+
+
+
+
    emptyOrTextOrCellOfText
+
    
+
+
    PURPOSE ^
    
+
    Validate [] OR text scalar (char row or string scalar) OR cell array of such text.
    
+
+
    SYNOPSIS ^
    
+
    function mustBeEmptyOrTextOrCellOfText(x) 
    
+
+
    DESCRIPTION ^
    
+
     Validate [] OR text scalar (char row or string scalar) OR cell array of such text.
    
+
+
+
    CROSS-REFERENCE INFORMATION ^
    
+This function calls:
+
      
+
    
+This function is called by:
+
      
+
    
+
+
+
+
+
    SOURCE CODE ^
    
+
    0001 function mustBeEmptyOrTextOrCellOfText(x)
+0002 % Validate [] OR text scalar (char row or string scalar) OR cell array of such text.
+0003 
+0004 % Allow explicit empty
+0005 if isempty(x)
+0006     return;
+0007 end
+0008 
+0009 % Allow a single text scalar
+0010 if ischar(x)
+0011     if ~isrow(x)
+0012         error('Text must be a character row vector, a string scalar, or empty.');
+0013     end
+0014     return;
+0015 elseif isstring(x)
+0016     if ~isscalar(x)
+0017         error('Text must be a string scalar, a character row vector, or empty.');
+0018     end
+0019     return;
+0020 end
+0021 
+0022 % Allow a cell array of text
+0023 if iscell(x)
+0024     for k = 1:numel(x)
+0025         v = x{k};
+0026         if ischar(v)
+0027             if ~isrow(v)   % '' is 1x0 char, which is a row -> allowed
+0028                 error('Each cell element must be a character row vector or a string scalar.');
+0029             end
+0030         elseif isstring(v)
+0031             if ~isscalar(v)  % "" is a scalar string -> allowed
+0032                 error('Each cell element must be a string scalar or a character row vector.');
+0033             end
+0034         else
+0035             error('Each cell element must be text (string scalar or character row vector).');
+0036         end
+0037     end
+0038     return;
+0039 end
+0040 
+0041 % Anything else is invalid
+0042 error('Value must be empty, a text value, or a cell array of text values.');
+0043 end
    
+
    Generated by m2html © 2005
    
+
+
\ No newline at end of file
diff --git a/doc/utils/emptyOrTextScalar.html b/doc/utils/emptyOrTextScalar.html
new file mode 100644
index 00000000..7d53acd3
--- /dev/null
+++ b/doc/utils/emptyOrTextScalar.html
@@ -0,0 +1,59 @@
+
+
+
+  Description of emptyOrTextScalar
+  
+  
+  
+  
+  
+  
+
+
+
+
    Home >  utils > emptyOrTextScalar.m
    
+
+
+
+
    emptyOrTextScalar
+
    
+
+
    PURPOSE ^
    
+
    
+
+
    SYNOPSIS ^
    
+
    function mustBeEmptyOrTextScalar(x) 
    
+
+
    DESCRIPTION ^
    
+
    
+
+
+
    CROSS-REFERENCE INFORMATION ^
    
+This function calls:
+
      
+
    
+This function is called by:
+
      
+
    
+
+
+
+
+
    SOURCE CODE ^
    
+
    0001 function mustBeEmptyOrTextScalar(x)
+0002     if isempty(x)
+0003         return;
+0004     end
+0005     if isstring(x)
+0006         if ~isscalar(x), error('Value must be a string scalar or empty.'); end
+0007     elseif ischar(x)
+0008         if ~isrow(x), error('Char input must be a row vector or empty.'); end
+0009     else
+0010         error('Value must be a string/char scalar or empty.');
+0011     end
+0012 end
    
+
    Generated by m2html © 2005
    
+
+
\ No newline at end of file
diff --git a/doc/utils/index.html b/doc/utils/index.html
new file mode 100644
index 00000000..f71a3bf0
--- /dev/null
+++ b/doc/utils/index.html
@@ -0,0 +1,29 @@
+
+
+
+  Index for Directory utils
+  
+  
+  
+  
+  
+  
+
+
+
+
+
    < Master indexIndex for utils >
    
+
+
    Index for utils
    
+
+
    Matlab files in this directory:
    
+
+
     emptyOrLogicalScalar 
     emptyOrTextOrCellOfTextValidate [] OR text scalar (char row or string scalar) OR cell array of such text. 
     emptyOrTextScalar 
    
+
+
+
+
+
    Generated by m2html © 2005
    
+
+
\ No newline at end of file
diff --git a/testing/unit_tests/importExportTests.m b/testing/unit_tests/importExportTests.m
index 4db6eb9f..ff670f57 100755
--- a/testing/unit_tests/importExportTests.m
+++ b/testing/unit_tests/importExportTests.m
@@ -36,7 +36,7 @@ function testYAMLimport(testCase)
 function testExcelExport(testCase)
 sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
 load(fullfile(sourceDir,'testing','unit_tests','test_data','ecoli_textbook.mat'), 'model');
-exportToExcelFormat(model,fullfile(sourceDir,'testing','unit_tests','test_data','_test.xlsx'));
+evalc('exportToExcelFormat(model,fullfile(sourceDir,''testing'',''unit_tests'',''test_data'',''_test.xlsx''))');
 %File will not be exactly equal as it contains the current date and time,
 %so md5 or similar would not work. Just check whether file is reasonably
 %sized.
diff --git a/utils/emptyOrLogicalScalar.m b/utils/emptyOrLogicalScalar.m
new file mode 100644
index 00000000..1583901d
--- /dev/null
+++ b/utils/emptyOrLogicalScalar.m
@@ -0,0 +1,6 @@
+function emptyOrLogical(x)
+    if isempty(x)
+        return;
+    end
+    validateattributes(x, {'logical'}, {'scalar'});
+end
diff --git a/utils/emptyOrTextOrCellOfText.m b/utils/emptyOrTextOrCellOfText.m
new file mode 100644
index 00000000..7dc95d5f
--- /dev/null
+++ b/utils/emptyOrTextOrCellOfText.m
@@ -0,0 +1,43 @@
+function mustBeEmptyOrTextOrCellOfText(x)
+% Validate [] OR text scalar (char row or string scalar) OR cell array of such text.
+
+% Allow explicit empty
+if isempty(x)
+    return;
+end
+
+% Allow a single text scalar
+if ischar(x)
+    if ~isrow(x)
+        error('Text must be a character row vector, a string scalar, or empty.');
+    end
+    return;
+elseif isstring(x)
+    if ~isscalar(x)
+        error('Text must be a string scalar, a character row vector, or empty.');
+    end
+    return;
+end
+
+% Allow a cell array of text
+if iscell(x)
+    for k = 1:numel(x)
+        v = x{k};
+        if ischar(v)
+            if ~isrow(v)   % '' is 1x0 char, which is a row -> allowed
+                error('Each cell element must be a character row vector or a string scalar.');
+            end
+        elseif isstring(v)
+            if ~isscalar(v)  % "" is a scalar string -> allowed
+                error('Each cell element must be a string scalar or a character row vector.');
+            end
+        else
+            error('Each cell element must be text (string scalar or character row vector).');
+        end
+    end
+    return;
+end
+
+% Anything else is invalid
+error('Value must be empty, a text value, or a cell array of text values.');
+end
diff --git a/utils/emptyOrTextScalar.m b/utils/emptyOrTextScalar.m
new file mode 100644
index 00000000..3c61cdd7
--- /dev/null
+++ b/utils/emptyOrTextScalar.m
@@ -0,0 +1,12 @@
+function mustBeEmptyOrTextScalar(x)
+    if isempty(x)
+        return;
+    end
+    if isstring(x)
+        if ~isscalar(x), error('Value must be a string scalar or empty.'); end
+    elseif ischar(x)
+        if ~isrow(x), error('Char input must be a row vector or empty.'); end
+    else
+        error('Value must be a string/char scalar or empty.');
+    end
+end

From 722e93cffcede3b3ca89cfa05946f8fe021caec8 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven 
Date: Thu, 15 Jan 2026 22:49:38 +0100
Subject: [PATCH 7/9] fix: addRxns make subSystems nested cell arrays

---
 core/addRxns.m        |  14 +-
 doc/core/addRxns.html | 628 +++++++++++++++++++++---------------------
 2 files changed, 325 insertions(+), 317 deletions(-)

diff --git a/core/addRxns.m b/core/addRxns.m
index f23222e1..8efcd5e2 100755
--- a/core/addRxns.m
+++ b/core/addRxns.m
@@ -350,15 +350,19 @@
 end
 
 if isfield(rxnsToAdd,'subSystems')
+    if ischar(rxnsToAdd.subSystems)
+        rxnsToAdd.subSystems = {{rxnsToAdd.subSystems}};
+    else
+        for i=1:numel(rxnsToAdd.subSystems)
+            if ischar(rxnsToAdd.subSystems{i})
+                rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
+            end
+        end
+    end
     if numel(rxnsToAdd.subSystems)~=nRxns
         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
         dispEM(EM);
     end
-    for i=1:numel(rxnsToAdd.subSystems)
-        if ischar(rxnsToAdd.subSystems{i})
-            rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
-        end
-    end
     %Fill with standard if it doesn't exist
     if ~isfield(newModel,'subSystems')
         newModel.subSystems=celllargefiller;
diff --git a/doc/core/addRxns.html b/doc/core/addRxns.html
index ffdcf6e6..e4ebe1c7 100644
--- a/doc/core/addRxns.html
+++ b/doc/core/addRxns.html
@@ -494,330 +494,334 @@ 
    SOURCE CODE ^end
 0351 
 0352 if isfield(rxnsToAdd,'subSystems')
-0353     if numel(rxnsToAdd.subSystems)~=nRxns
-0354         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
-0355         dispEM(EM);
-0356     end
-0357     for i=1:numel(rxnsToAdd.subSystems)
-0358         if ischar(rxnsToAdd.subSystems{i})
-0359             rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
+0353     if ischar(rxnsToAdd.subSystems)
+0354         rxnsToAdd.subSystems = {{rxnsToAdd.subSystems}};
+0355     else
+0356         for i=1:numel(rxnsToAdd.subSystems)
+0357             if ischar(rxnsToAdd.subSystems{i})
+0358                 rxnsToAdd.subSystems{i}=rxnsToAdd.subSystems(i);
+0359             end
 0360         end
 0361     end
-0362     %Fill with standard if it doesn't exist
-0363     if ~isfield(newModel,'subSystems')
-0364         newModel.subSystems=celllargefiller;
+0362     if numel(rxnsToAdd.subSystems)~=nRxns
+0363         EM='rxnsToAdd.subSystems must have the same number of elements as rxnsToAdd.rxns';
+0364         dispEM(EM);
 0365     end
-0366     newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
-0367 else
-0368     %Fill with standard if it doesn't exist
-0369     if isfield(newModel,'subSystems')
-0370         newModel.subSystems=[newModel.subSystems;cellfiller];
-0371     end
-0372 end
-0373 if isfield(rxnsToAdd,'rxnMiriams')
-0374     if numel(rxnsToAdd.rxnMiriams)~=nRxns
-0375         EM='rxnsToAdd.rxnMiriams must have the same number of elements as rxnsToAdd.rxns';
-0376         dispEM(EM);
-0377     end
-0378     %Fill with standard if it doesn't exist
-0379     if ~isfield(newModel,'rxnMiriams')
-0380         newModel.rxnMiriams=cell(nOldRxns,1);
+0366     %Fill with standard if it doesn't exist
+0367     if ~isfield(newModel,'subSystems')
+0368         newModel.subSystems=celllargefiller;
+0369     end
+0370     newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)];
+0371 else
+0372     %Fill with standard if it doesn't exist
+0373     if isfield(newModel,'subSystems')
+0374         newModel.subSystems=[newModel.subSystems;cellfiller];
+0375     end
+0376 end
+0377 if isfield(rxnsToAdd,'rxnMiriams')
+0378     if numel(rxnsToAdd.rxnMiriams)~=nRxns
+0379         EM='rxnsToAdd.rxnMiriams must have the same number of elements as rxnsToAdd.rxns';
+0380         dispEM(EM);
 0381     end
-0382     newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
-0383 else
-0384     %Fill with standard if it doesn't exist
-0385     if isfield(newModel,'rxnMiriams')
-0386         newModel.rxnMiriams=[newModel.rxnMiriams;cell(nRxns,1)];
-0387     end
-0388 end
-0389 if isfield(rxnsToAdd,'rxnComps')
-0390     if numel(rxnsToAdd.rxnComps)~=nRxns
-0391         EM='rxnsToAdd.rxnComps must have the same number of elements as rxnsToAdd.rxns';
-0392         dispEM(EM);
-0393     end
-0394     %Fill with standard if it doesn't exist
-0395     if ~isfield(newModel,'rxnComps')
-0396         newModel.rxnComps=ones(nOldRxns,1);
-0397         EM='Adding reactions with compartment information to a model without such information. All existing reactions will be assigned to the first compartment';
-0398         dispEM(EM,false);
-0399     end
-0400     newModel.rxnComps=[newModel.rxnComps;rxnsToAdd.rxnComps(:)];
-0401 else
-0402     %Fill with standard if it doesn't exist
-0403     if isfield(newModel,'rxnComps')
-0404         newModel.rxnComps=[newModel.rxnComps;ones(nRxns,1)];
-0405         %fprintf('NOTE: The added reactions will be assigned to the first
-0406         %compartment\n');
-0407     end
-0408 end
-0409 if isfield(rxnsToAdd,'grRules')
-0410     rxnsToAdd.grRules=convertCharArray(rxnsToAdd.grRules);
-0411     if numel(rxnsToAdd.grRules)~=nRxns
-0412         EM='rxnsToAdd.grRules must have the same number of elements as rxnsToAdd.rxns';
-0413         dispEM(EM);
-0414     end
-0415     %Fill with standard if it doesn't exist
-0416     if ~isfield(newModel,'grRules')
-0417         newModel.grRules=largeFiller;
+0382     %Fill with standard if it doesn't exist
+0383     if ~isfield(newModel,'rxnMiriams')
+0384         newModel.rxnMiriams=cell(nOldRxns,1);
+0385     end
+0386     newModel.rxnMiriams=[newModel.rxnMiriams;rxnsToAdd.rxnMiriams(:)];
+0387 else
+0388     %Fill with standard if it doesn't exist
+0389     if isfield(newModel,'rxnMiriams')
+0390         newModel.rxnMiriams=[newModel.rxnMiriams;cell(nRxns,1)];
+0391     end
+0392 end
+0393 if isfield(rxnsToAdd,'rxnComps')
+0394     if numel(rxnsToAdd.rxnComps)~=nRxns
+0395         EM='rxnsToAdd.rxnComps must have the same number of elements as rxnsToAdd.rxns';
+0396         dispEM(EM);
+0397     end
+0398     %Fill with standard if it doesn't exist
+0399     if ~isfield(newModel,'rxnComps')
+0400         newModel.rxnComps=ones(nOldRxns,1);
+0401         EM='Adding reactions with compartment information to a model without such information. All existing reactions will be assigned to the first compartment';
+0402         dispEM(EM,false);
+0403     end
+0404     newModel.rxnComps=[newModel.rxnComps;rxnsToAdd.rxnComps(:)];
+0405 else
+0406     %Fill with standard if it doesn't exist
+0407     if isfield(newModel,'rxnComps')
+0408         newModel.rxnComps=[newModel.rxnComps;ones(nRxns,1)];
+0409         %fprintf('NOTE: The added reactions will be assigned to the first
+0410         %compartment\n');
+0411     end
+0412 end
+0413 if isfield(rxnsToAdd,'grRules')
+0414     rxnsToAdd.grRules=convertCharArray(rxnsToAdd.grRules);
+0415     if numel(rxnsToAdd.grRules)~=nRxns
+0416         EM='rxnsToAdd.grRules must have the same number of elements as rxnsToAdd.rxns';
+0417         dispEM(EM);
 0418     end
-0419     newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
-0420 else
-0421     %Fill with standard if it doesn't exist
-0422     if isfield(newModel,'grRules')
-0423         newModel.grRules=[newModel.grRules;filler];
-0424     end
-0425 end
-0426 
-0427 if isfield(rxnsToAdd,'rxnFrom')
-0428     rxnsToAdd.rxnFrom=convertCharArray(rxnsToAdd.rxnFrom);
-0429     if numel(rxnsToAdd.rxnFrom)~=nRxns
-0430         EM='rxnsToAdd.rxnFrom must have the same number of elements as rxnsToAdd.rxns';
-0431         dispEM(EM);
-0432     end
-0433     %Fill with standard if it doesn't exist
-0434     if ~isfield(newModel,'rxnFrom')
-0435         newModel.rxnFrom=largeFiller;
+0419     %Fill with standard if it doesn't exist
+0420     if ~isfield(newModel,'grRules')
+0421         newModel.grRules=largeFiller;
+0422     end
+0423     newModel.grRules=[newModel.grRules;rxnsToAdd.grRules(:)];
+0424 else
+0425     %Fill with standard if it doesn't exist
+0426     if isfield(newModel,'grRules')
+0427         newModel.grRules=[newModel.grRules;filler];
+0428     end
+0429 end
+0430 
+0431 if isfield(rxnsToAdd,'rxnFrom')
+0432     rxnsToAdd.rxnFrom=convertCharArray(rxnsToAdd.rxnFrom);
+0433     if numel(rxnsToAdd.rxnFrom)~=nRxns
+0434         EM='rxnsToAdd.rxnFrom must have the same number of elements as rxnsToAdd.rxns';
+0435         dispEM(EM);
 0436     end
-0437     newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
-0438 else
-0439     %Fill with standard if it doesn't exist
-0440     if isfield(newModel,'rxnFrom')
-0441         newModel.rxnFrom=[newModel.rxnFrom;filler];
-0442     end
-0443 end
-0444 
-0445 if isfield(rxnsToAdd,'rxnNotes')
-0446     rxnsToAdd.rxnNotes=convertCharArray(rxnsToAdd.rxnNotes);
-0447     if numel(rxnsToAdd.rxnNotes)~=nRxns
-0448         EM='rxnsToAdd.rxnNotes must have the same number of elements as rxnsToAdd.rxns';
-0449         dispEM(EM);
-0450     end
-0451     %Fill with standard if it doesn't exist
-0452     if ~isfield(newModel,'rxnNotes')
-0453         newModel.rxnNotes=largeFiller;
+0437     %Fill with standard if it doesn't exist
+0438     if ~isfield(newModel,'rxnFrom')
+0439         newModel.rxnFrom=largeFiller;
+0440     end
+0441     newModel.rxnFrom=[newModel.rxnFrom;rxnsToAdd.rxnFrom(:)];
+0442 else
+0443     %Fill with standard if it doesn't exist
+0444     if isfield(newModel,'rxnFrom')
+0445         newModel.rxnFrom=[newModel.rxnFrom;filler];
+0446     end
+0447 end
+0448 
+0449 if isfield(rxnsToAdd,'rxnNotes')
+0450     rxnsToAdd.rxnNotes=convertCharArray(rxnsToAdd.rxnNotes);
+0451     if numel(rxnsToAdd.rxnNotes)~=nRxns
+0452         EM='rxnsToAdd.rxnNotes must have the same number of elements as rxnsToAdd.rxns';
+0453         dispEM(EM);
 0454     end
-0455     newModel.rxnNotes=[newModel.rxnNotes;rxnsToAdd.rxnNotes(:)];
-0456 else
-0457     %Fill with standard if it doesn't exist
-0458     if isfield(newModel,'rxnNotes')
-0459         newModel.rxnNotes=[newModel.rxnNotes;filler];
-0460     end
-0461 end
-0462 
-0463 if isfield(rxnsToAdd,'rxnReferences')
-0464     rxnsToAdd.rxnReferences=convertCharArray(rxnsToAdd.rxnReferences);
-0465     if numel(rxnsToAdd.rxnReferences)~=nRxns
-0466         EM='rxnsToAdd.rxnReferences must have the same number of elements as rxnsToAdd.rxns';
-0467         dispEM(EM);
-0468     end
-0469     %Fill with standard if it doesn't exist
-0470     if ~isfield(newModel,'rxnReferences')
-0471         newModel.rxnReferences=largeFiller;
+0455     %Fill with standard if it doesn't exist
+0456     if ~isfield(newModel,'rxnNotes')
+0457         newModel.rxnNotes=largeFiller;
+0458     end
+0459     newModel.rxnNotes=[newModel.rxnNotes;rxnsToAdd.rxnNotes(:)];
+0460 else
+0461     %Fill with standard if it doesn't exist
+0462     if isfield(newModel,'rxnNotes')
+0463         newModel.rxnNotes=[newModel.rxnNotes;filler];
+0464     end
+0465 end
+0466 
+0467 if isfield(rxnsToAdd,'rxnReferences')
+0468     rxnsToAdd.rxnReferences=convertCharArray(rxnsToAdd.rxnReferences);
+0469     if numel(rxnsToAdd.rxnReferences)~=nRxns
+0470         EM='rxnsToAdd.rxnReferences must have the same number of elements as rxnsToAdd.rxns';
+0471         dispEM(EM);
 0472     end
-0473     newModel.rxnReferences=[newModel.rxnReferences;rxnsToAdd.rxnReferences(:)];
-0474 else
-0475     %Fill with standard if it doesn't exist
-0476     if isfield(newModel,'rxnReferences')
-0477         newModel.rxnReferences=[newModel.rxnReferences;filler];
-0478     end
-0479 end
-0480 
-0481 if isfield(rxnsToAdd,'pwys')
-0482     rxnsToAdd.pwys=convertCharArray(rxnsToAdd.pwys);
-0483     if numel(rxnsToAdd.pwys)~=nRxns
-0484         EM='rxnsToAdd.pwys must have the same number of elements as rxnsToAdd.rxns';
-0485         dispEM(EM);
-0486     end
-0487     %Fill with standard if it doesn't exist
-0488     if ~isfield(newModel,'pwys')
-0489         newModel.pwys=largeFiller;
+0473     %Fill with standard if it doesn't exist
+0474     if ~isfield(newModel,'rxnReferences')
+0475         newModel.rxnReferences=largeFiller;
+0476     end
+0477     newModel.rxnReferences=[newModel.rxnReferences;rxnsToAdd.rxnReferences(:)];
+0478 else
+0479     %Fill with standard if it doesn't exist
+0480     if isfield(newModel,'rxnReferences')
+0481         newModel.rxnReferences=[newModel.rxnReferences;filler];
+0482     end
+0483 end
+0484 
+0485 if isfield(rxnsToAdd,'pwys')
+0486     rxnsToAdd.pwys=convertCharArray(rxnsToAdd.pwys);
+0487     if numel(rxnsToAdd.pwys)~=nRxns
+0488         EM='rxnsToAdd.pwys must have the same number of elements as rxnsToAdd.rxns';
+0489         dispEM(EM);
 0490     end
-0491     newModel.pwys=[newModel.pwys;rxnsToAdd.pwys(:)];
-0492 else
-0493     %Fill with standard if it doesn't exist
-0494     if isfield(newModel,'pwys')
-0495         newModel.pwys=[newModel.pwys;filler];
-0496     end
-0497 end
-0498 
-0499 if isfield(rxnsToAdd,'rxnConfidenceScores')
-0500     if numel(rxnsToAdd.rxnConfidenceScores)~=nRxns
-0501         EM='rxnsToAdd.rxnConfidenceScores must have the same number of elements as rxnsToAdd.rxns';
-0502         dispEM(EM);
-0503     end
-0504     %Fill with standard if it doesn't exist
-0505     if ~isfield(newModel,'rxnConfidenceScores')
-0506         newModel.rxnConfidenceScores=NaN(nOldRxns,1);
+0491     %Fill with standard if it doesn't exist
+0492     if ~isfield(newModel,'pwys')
+0493         newModel.pwys=largeFiller;
+0494     end
+0495     newModel.pwys=[newModel.pwys;rxnsToAdd.pwys(:)];
+0496 else
+0497     %Fill with standard if it doesn't exist
+0498     if isfield(newModel,'pwys')
+0499         newModel.pwys=[newModel.pwys;filler];
+0500     end
+0501 end
+0502 
+0503 if isfield(rxnsToAdd,'rxnConfidenceScores')
+0504     if numel(rxnsToAdd.rxnConfidenceScores)~=nRxns
+0505         EM='rxnsToAdd.rxnConfidenceScores must have the same number of elements as rxnsToAdd.rxns';
+0506         dispEM(EM);
 0507     end
-0508     newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;rxnsToAdd.rxnConfidenceScores(:)];
-0509 else
-0510     %Fill with standard if it doesn't exist
-0511     if isfield(newModel,'rxnConfidenceScores')
-0512         newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;NaN(nRxns,1)];
-0513     end
-0514 end
-0515 
-0516 if isfield(rxnsToAdd,'rxnDeltaG')
-0517     if numel(rxnsToAdd.rxnDeltaG)~=nRxns
-0518         EM='rxnsToAdd.rxnDeltaG must have the same number of elements as rxnsToAdd.rxns';
-0519         dispEM(EM);
-0520     end
-0521     %Fill with standard if it doesn't exist
-0522     if ~isfield(newModel,'rxnDeltaG')
-0523         newModel.rxnDeltaG=NaN(nOldRxns,1);
+0508     %Fill with standard if it doesn't exist
+0509     if ~isfield(newModel,'rxnConfidenceScores')
+0510         newModel.rxnConfidenceScores=NaN(nOldRxns,1);
+0511     end
+0512     newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;rxnsToAdd.rxnConfidenceScores(:)];
+0513 else
+0514     %Fill with standard if it doesn't exist
+0515     if isfield(newModel,'rxnConfidenceScores')
+0516         newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;NaN(nRxns,1)];
+0517     end
+0518 end
+0519 
+0520 if isfield(rxnsToAdd,'rxnDeltaG')
+0521     if numel(rxnsToAdd.rxnDeltaG)~=nRxns
+0522         EM='rxnsToAdd.rxnDeltaG must have the same number of elements as rxnsToAdd.rxns';
+0523         dispEM(EM);
 0524     end
-0525     newModel.rxnDeltaG=[newModel.rxnDeltaG;rxnsToAdd.rxnDeltaG(:)];
-0526 else
-0527     %Fill with standard if it doesn't exist
-0528     if isfield(newModel,'rxnDeltaG')
-0529         newModel.rxnDeltaG=[newModel.rxnDeltaG;NaN(nRxns,1)];
-0530     end
-0531 end
-0532 
-0533 
-0534 %***Start parsing the equations and adding the info to the S matrix The
-0535 %mets are matched to model.mets
-0536 if eqnType==1
-0537     [I, J]=ismember(mets,model.mets);
-0538     if ~all(I)
-0539         if allowNewMets==true || ischar(allowNewMets)
-0540             %Add the new mets
-0541             metsToAdd.mets=mets(~I);
-0542             metsToAdd.metNames=metsToAdd.mets;
-0543             metsToAdd.compartments=compartment;
-0544             if ischar(allowNewMets)
-0545                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
-0546             else
-0547                 newModel=addMets(newModel,metsToAdd,true);
-0548             end
-0549         else
-0550             EM='One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function. Are you sure that eqnType=1?';
-0551             dispEM(EM);
-0552         end
-0553     end
-0554     %Calculate the indexes of the metabolites and add the info
-0555     metIndexes=J;
-0556     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
-0557 end
-0558 
-0559 %Do some stuff that is the same for eqnType=2 and eqnType=3
-0560 if eqnType==2 || eqnType==3
-0561     %For later..
-0562     t2=strcat(model.metNames,'***',model.comps(model.metComps));
-0563 end
-0564 
-0565 %The mets are matched to model.metNames and assigned to "compartment"
-0566 if eqnType==2
-0567     %%Check that the metabolite names aren't present in the same
-0568     %%compartment.
-0569     %Not the neatest way maybe..
-0570     t1=strcat(mets,'***',compartment);
-0571     [I, J]=ismember(t1,t2);
-0572     
-0573     if ~all(I)
-0574         if allowNewMets==true || ischar(allowNewMets)
-0575             %Add the new mets
-0576             metsToAdd.metNames=mets(~I);
-0577             metsToAdd.compartments=compartment;
-0578             if ischar(allowNewMets)
-0579                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
-0580             else
-0581                 newModel=addMets(newModel,metsToAdd,true);
-0582             end
-0583         else
-0584             EM='One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function';
-0585             dispEM(EM);
-0586         end
-0587     end
-0588     
-0589     %Calculate the indexes of the metabolites
-0590     metIndexes=J;
-0591     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
-0592 end
-0593 
-0594 %The equations are on the form metNames[compName]
-0595 if eqnType==3
-0596     %Parse the metabolite names
-0597     metNames=cell(numel(mets),1);
-0598     compartments=metNames;
-0599     for i=1:numel(mets)
-0600         starts=max(strfind(mets{i},'['));
-0601         ends=max(strfind(mets{i},']'));
-0602         
-0603         %Check that the formatting is correct
-0604         if isempty(starts) || isempty(ends) || ends<numel(mets{i})
-0605             EM=['The metabolite ' mets{i} ' is not correctly formatted for eqnType=3'];
-0606             dispEM(EM);
-0607         end
-0608         
-0609         %Check that the compartment is correct
-0610         compartments{i}=mets{i}(starts+1:ends-1);
-0611         I=ismember(compartments(i),newModel.comps);
-0612         if ~I
-0613             EM=['The metabolite ' mets{i} ' has a compartment that is not in model.comps'];
-0614             dispEM(EM);
-0615         end
-0616         metNames{i}=mets{i}(1:starts-1);
-0617     end
-0618     
-0619     %Check if the metabolite exists already
-0620     t1=strcat(metNames,'***',compartments);
-0621     [I, J]=ismember(t1,t2);
+0525     %Fill with standard if it doesn't exist
+0526     if ~isfield(newModel,'rxnDeltaG')
+0527         newModel.rxnDeltaG=NaN(nOldRxns,1);
+0528     end
+0529     newModel.rxnDeltaG=[newModel.rxnDeltaG;rxnsToAdd.rxnDeltaG(:)];
+0530 else
+0531     %Fill with standard if it doesn't exist
+0532     if isfield(newModel,'rxnDeltaG')
+0533         newModel.rxnDeltaG=[newModel.rxnDeltaG;NaN(nRxns,1)];
+0534     end
+0535 end
+0536 
+0537 
+0538 %***Start parsing the equations and adding the info to the S matrix The
+0539 %mets are matched to model.mets
+0540 if eqnType==1
+0541     [I, J]=ismember(mets,model.mets);
+0542     if ~all(I)
+0543         if allowNewMets==true || ischar(allowNewMets)
+0544             %Add the new mets
+0545             metsToAdd.mets=mets(~I);
+0546             metsToAdd.metNames=metsToAdd.mets;
+0547             metsToAdd.compartments=compartment;
+0548             if ischar(allowNewMets)
+0549                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
+0550             else
+0551                 newModel=addMets(newModel,metsToAdd,true);
+0552             end
+0553         else
+0554             EM='One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function. Are you sure that eqnType=1?';
+0555             dispEM(EM);
+0556         end
+0557     end
+0558     %Calculate the indexes of the metabolites and add the info
+0559     metIndexes=J;
+0560     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
+0561 end
+0562 
+0563 %Do some stuff that is the same for eqnType=2 and eqnType=3
+0564 if eqnType==2 || eqnType==3
+0565     %For later..
+0566     t2=strcat(model.metNames,'***',model.comps(model.metComps));
+0567 end
+0568 
+0569 %The mets are matched to model.metNames and assigned to "compartment"
+0570 if eqnType==2
+0571     %%Check that the metabolite names aren't present in the same
+0572     %%compartment.
+0573     %Not the neatest way maybe..
+0574     t1=strcat(mets,'***',compartment);
+0575     [I, J]=ismember(t1,t2);
+0576     
+0577     if ~all(I)
+0578         if allowNewMets==true || ischar(allowNewMets)
+0579             %Add the new mets
+0580             metsToAdd.metNames=mets(~I);
+0581             metsToAdd.compartments=compartment;
+0582             if ischar(allowNewMets)
+0583                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
+0584             else
+0585                 newModel=addMets(newModel,metsToAdd,true);
+0586             end
+0587         else
+0588             EM='One or more equations contain metabolites that are not in model.mets. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function';
+0589             dispEM(EM);
+0590         end
+0591     end
+0592     
+0593     %Calculate the indexes of the metabolites
+0594     metIndexes=J;
+0595     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
+0596 end
+0597 
+0598 %The equations are on the form metNames[compName]
+0599 if eqnType==3
+0600     %Parse the metabolite names
+0601     metNames=cell(numel(mets),1);
+0602     compartments=metNames;
+0603     for i=1:numel(mets)
+0604         starts=max(strfind(mets{i},'['));
+0605         ends=max(strfind(mets{i},']'));
+0606         
+0607         %Check that the formatting is correct
+0608         if isempty(starts) || isempty(ends) || ends<numel(mets{i})
+0609             EM=['The metabolite ' mets{i} ' is not correctly formatted for eqnType=3'];
+0610             dispEM(EM);
+0611         end
+0612         
+0613         %Check that the compartment is correct
+0614         compartments{i}=mets{i}(starts+1:ends-1);
+0615         I=ismember(compartments(i),newModel.comps);
+0616         if ~I
+0617             EM=['The metabolite ' mets{i} ' has a compartment that is not in model.comps'];
+0618             dispEM(EM);
+0619         end
+0620         metNames{i}=mets{i}(1:starts-1);
+0621     end
 0622     
-0623     if ~all(I)
-0624         if allowNewMets==true | ischar(allowNewMets)
-0625             %Add the new mets
-0626             metsToAdd.metNames=metNames(~I);
-0627             metsToAdd.compartments=compartments(~I);
-0628             if ischar(allowNewMets)
-0629                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
-0630             else
-0631                 newModel=addMets(newModel,metsToAdd,true);
-0632             end
-0633         else
-0634             EM='One or more equations contain metabolites that are not in model.metNames. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function';
-0635             dispEM(EM);
-0636         end
-0637     end
-0638     
-0639     %Calculate the indexes of the metabolites
-0640     metIndexes=J;
-0641     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
-0642 end
-0643 
-0644 %Add the info to the stoichiometric matrix
-0645 newModel.S=[newModel.S sparse(size(newModel.S,1),nRxns)];
-0646 
-0647 for i=1:nRxns
-0648     newModel.S(metIndexes,nOldRxns+i)=S(:,i);
-0649     %Parse the grRules and check whether all genes in grRules appear in
-0650     %genes
-0651     if isfield(newModel,'grRules')
-0652         rule=newModel.grRules{nOldRxns+i};
-0653         rule=strrep(rule,'(','');
-0654         rule=strrep(rule,')','');
-0655         rule=strrep(rule,' or ',' ');
-0656         rule=strrep(rule,' and ',' ');
-0657         genes=regexp(rule,' ','split');
-0658         [I, J]=ismember(genes,newModel.genes);
-0659         if ~all(I) && any(rule)
-0660             EM=['Not all genes for reaction ' rxnsToAdd.rxns{i} ' were found in model.genes. If needed, add genes with addGenesRaven before calling this function, or set the ''allowNewGenes'' flag to true'];
-0661             dispEM(EM);
-0662         end
-0663     end
-0664 end
-0665 
-0666 %Make temporary minimal model structure with only new rxns, to parse to
-0667 %standardizeGrRules
-0668 newRxnsModel.genes=newModel.genes;
-0669 newRxnsModel.grRules=newModel.grRules(length(model.rxns)+1:end);
-0670 newRxnsModel.rxns=newModel.rxns(length(model.rxns)+1:end);
-0671 
-0672 %Fix grRules and reconstruct rxnGeneMat
-0673 [grRules,rxnGeneMat] = standardizeGrRules(newRxnsModel,true);
-0674 newModel.rxnGeneMat = [newModel.rxnGeneMat; rxnGeneMat];
-0675 newModel.grRules = [newModel.grRules(1:nOldRxns); grRules];
-0676 end
    
+0623     %Check if the metabolite exists already
+0624     t1=strcat(metNames,'***',compartments);
+0625     [I, J]=ismember(t1,t2);
+0626     
+0627     if ~all(I)
+0628         if allowNewMets==true | ischar(allowNewMets)
+0629             %Add the new mets
+0630             metsToAdd.metNames=metNames(~I);
+0631             metsToAdd.compartments=compartments(~I);
+0632             if ischar(allowNewMets)
+0633                 newModel=addMets(newModel,metsToAdd,true,allowNewMets);
+0634             else
+0635                 newModel=addMets(newModel,metsToAdd,true);
+0636             end
+0637         else
+0638             EM='One or more equations contain metabolites that are not in model.metNames. Set allowNewMets to true to allow this function to add metabolites or use addMets to add them before calling this function';
+0639             dispEM(EM);
+0640         end
+0641     end
+0642     
+0643     %Calculate the indexes of the metabolites
+0644     metIndexes=J;
+0645     metIndexes(~I)=numel(newModel.mets)-sum(~I)+1:numel(newModel.mets);
+0646 end
+0647 
+0648 %Add the info to the stoichiometric matrix
+0649 newModel.S=[newModel.S sparse(size(newModel.S,1),nRxns)];
+0650 
+0651 for i=1:nRxns
+0652     newModel.S(metIndexes,nOldRxns+i)=S(:,i);
+0653     %Parse the grRules and check whether all genes in grRules appear in
+0654     %genes
+0655     if isfield(newModel,'grRules')
+0656         rule=newModel.grRules{nOldRxns+i};
+0657         rule=strrep(rule,'(','');
+0658         rule=strrep(rule,')','');
+0659         rule=strrep(rule,' or ',' ');
+0660         rule=strrep(rule,' and ',' ');
+0661         genes=regexp(rule,' ','split');
+0662         [I, J]=ismember(genes,newModel.genes);
+0663         if ~all(I) && any(rule)
+0664             EM=['Not all genes for reaction ' rxnsToAdd.rxns{i} ' were found in model.genes. If needed, add genes with addGenesRaven before calling this function, or set the ''allowNewGenes'' flag to true'];
+0665             dispEM(EM);
+0666         end
+0667     end
+0668 end
+0669 
+0670 %Make temporary minimal model structure with only new rxns, to parse to
+0671 %standardizeGrRules
+0672 newRxnsModel.genes=newModel.genes;
+0673 newRxnsModel.grRules=newModel.grRules(length(model.rxns)+1:end);
+0674 newRxnsModel.rxns=newModel.rxns(length(model.rxns)+1:end);
+0675 
+0676 %Fix grRules and reconstruct rxnGeneMat
+0677 [grRules,rxnGeneMat] = standardizeGrRules(newRxnsModel,true);
+0678 newModel.rxnGeneMat = [newModel.rxnGeneMat; rxnGeneMat];
+0679 newModel.grRules = [newModel.grRules(1:nOldRxns); grRules];
+0680 end
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file

From 027f90af59babd4f9204cf1bb19b01a23204eb96 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven 
Date: Thu, 15 Jan 2026 22:49:56 +0100
Subject: [PATCH 8/9] refactor: exportModel speedup subSystem parsing

---
 doc/io/exportModel.html | 491 ++++++++++++++++++++--------------------
 io/exportModel.m        | 127 ++++++-----
 2 files changed, 318 insertions(+), 300 deletions(-)

diff --git a/doc/io/exportModel.html b/doc/io/exportModel.html
index 57f67266..6367bf1a 100644
--- a/doc/io/exportModel.html
+++ b/doc/io/exportModel.html
@@ -308,7 +308,7 @@ 
    SOURCE CODE ^if i<numel(model.comps)
 0246         modelSBML.compartment(i+1)=modelSBML.compartment(i);
 0247     end
-0248     
+0248 
 0249     if isfield(modelSBML.compartment,'metaid')
 0250         modelSBML.compartment(i).metaid=model.comps{i};
 0251     end
@@ -334,7 +334,7 @@ 
    SOURCE CODE ^if isfield(modelSBML.compartment, 'id')
 0272         modelSBML.compartment(i).id=model.comps{i};
 0273     end
-0274     
+0274 
 0275 end
 0276 
 0277 %Begin writing species
@@ -362,7 +362,7 @@ 
    SOURCE CODE ^if i<numel(model.mets)
 0300         modelSBML.species(i+1)=modelSBML.species(i);
 0301     end
-0302     
+0302 
 0303     if isfield(modelSBML.species,'metaid')
 0304         modelSBML.species(i).metaid=model.mets{i};
 0305     end
@@ -440,7 +440,7 @@ 
    SOURCE CODE ^if i<numel(model.genes)
 0378             modelSBML.fbc_geneProduct(i+1)=modelSBML.fbc_geneProduct(i);
 0379         end
-0380         
+0380 
 0381         if isfield(modelSBML.fbc_geneProduct,'metaid')
 0382             modelSBML.fbc_geneProduct(i).metaid=model.genes{i};
 0383         end
@@ -523,11 +523,11 @@ 
    SOURCE CODE ^if i<numel(model.rxns)
 0461         modelSBML.reaction(i+1)=modelSBML.reaction(i);
 0462     end
-0463     
+0463 
 0464     if isfield(modelSBML.reaction,'metaid')
 0465         modelSBML.reaction(i).metaid=model.rxns{i};
 0466     end
-0467     
+0467 
 0468     %Export notes information
 0469     if (~isnan(model.rxnConfidenceScores(i)) || ~isempty(model.rxnReferences{i}) || ~isempty(model.rxnNotes{i}))
 0470         modelSBML.reaction(i).notes='<notes><body xmlns="http://www.w3.org/1999/xhtml">';
@@ -542,7 +542,7 @@ 
    SOURCE CODE ^end
 0480         modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '</body></notes>'];
 0481     end
-0482     
+0482 
 0483     % Export SBO terms from rxnMiriams
 0484     if ~isempty(model.rxnMiriams{i})
 0485         [~,sbo_ind] = ismember('sbo',model.rxnMiriams{i}.name);
@@ -554,7 +554,7 @@ 
    SOURCE CODE ^end
 0493     end
-0494     
+0494 
 0495     %Export annotation information from rxnMiriams
 0496     if (~isempty(model.rxnMiriams{i}) && isfield(modelSBML.reaction(i),'annotation')) || ~isempty(model.eccodes{i})
 0497         modelSBML.reaction(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.rxns{i} '">'];
@@ -567,14 +567,14 @@ 
    SOURCE CODE ^end
 0505         modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation getMiriam(model.rxnMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
 0506     end
-0507     
+0507 
 0508     if isfield(modelSBML.reaction, 'name')
 0509         modelSBML.reaction(i).name=model.rxnNames{i};
 0510     end
 0511     if isfield(modelSBML.reaction, 'id')
 0512         modelSBML.reaction(i).id=model.rxns{i};
 0513     end
-0514     
+0514 
 0515     %Add the information about reactants and products
 0516     involvedMets=addReactantsProducts(model,modelSBML,i);
 0517     for j=1:numel(involvedMets.reactant)
@@ -616,241 +616,250 @@ 
    SOURCE CODE ^end
 0555 
-0556 %Prepare subSystems Code taken from COBRA functions getModelSubSystems,
-0557 %writeSBML, findRxnsFromSubSystem under GNU General Public License v3.0,
-0558 %license file in readme/GPL.MD. Code modified for RAVEN
-0559 if modelHasSubsystems
-0560     modelSBML.groups_group.groups_kind = 'partonomy';
-0561     modelSBML.groups_group.sboTerm = 633;
-0562     tmpStruct=modelSBML.groups_group;
-0563 
-0564     rxns=model.rxns;
-0565     if ~any(cellfun(@iscell,model.subSystems))
-0566         if ~any(~cellfun(@isempty,model.subSystems))
-0567             subSystems = {};
-0568         else
-0569             subSystems = setdiff(model.subSystems,'');
-0570         end
-0571     else
-0572         orderedSubs = cellfun(@(x) columnVector(x),model.subSystems,'UniformOUtput',false);
-0573         subSystems = setdiff(vertcat(orderedSubs{:}),'');
-0574     end
-0575     if isempty(subSystems)
-0576         subSystems = {};
-0577     end
-0578     if ~isempty(subSystems)
-0579         %Build the groups for the group package
-0580         groupIDs = strcat('group',cellfun(@num2str, num2cell(1:length(subSystems)),'UniformOutput',false));
-0581         for i = 1:length(subSystems)
-0582             cgroup = tmpStruct;
-0583             if ~any(cellfun(@iscell,model.subSystems))
-0584                 present = ismember(model.subSystems,subSystems{i});
-0585             else
-0586                 present = cellfun(@(x) any(ismember(x,subSystems{i})),model.subSystems);
-0587             end
-0588             groupMembers = rxns(present);
-0589             for j = 1:numel(groupMembers)
-0590                 cMember = tmpStruct.groups_member;
-0591                 cMember.groups_idRef = groupMembers{j};
-0592                 if j == 1
-0593                     cgroup.groups_member = cMember;
-0594                 else
-0595                     cgroup.groups_member(j) = cMember;
-0596                 end
-0597             end
-0598             cgroup.groups_id = groupIDs{i};
-0599             cgroup.groups_name = subSystems{i};
-0600             if i == 1
-0601                 modelSBML.groups_group = cgroup;
-0602             else
-0603                 modelSBML.groups_group(i) = cgroup;
-0604             end
-0605         end
-0606     end
-0607 end
-0608 
-0609 %Prepare fbc_objective subfield
+0556 %Prepare subSystems
+0557 if modelHasSubsystems
+0558     modelSBML.groups_group.groups_kind = 'partonomy';
+0559     modelSBML.groups_group.sboTerm     = 633;
+0560     grpTemplate = modelSBML.groups_group;
+0561 
+0562     % === 1) Normalize: make every entry a cell array of chars (vectorized) ===
+0563     isChar = cellfun(@ischar, model.subSystems);
+0564     model.subSystems(isChar)  = cellfun(@(s){s}, model.subSystems(isChar),  'UniformOutput', false);
+0565     model.subSystems(cellfun(@isempty, model.subSystems)) = {{}};
+0566 
+0567     % If some entries contain string scalars by mistake, coerce them:
+0568     model.subSystems = cellfun(@(c) cellfun(@char, c, 'UniformOutput', false), model.subSystems, 'UniformOutput', false);
+0569 
+0570     % === 2) Flatten once: names and their reaction indices (vectorized) ===
+0571     flatNames = [model.subSystems{:}];                    % 1×M cellstr of all subsystem labels
+0572     if isempty(flatNames)
+0573         % Nothing to do: no subsystems present
+0574         return
+0575     end
+0576 
+0577     counts   = cellfun(@numel, model.subSystems);         % reactions -> how many subsystems
+0578     % For each reaction r, repeat r exactly counts(r) times
+0579     flatIdx  = arrayfun(@(r,c) repmat(r, c, 1), (1:numel(model.subSystems)).', counts, 'UniformOutput', false);
+0580     flatIdx  = vertcat(flatIdx{:});           % M×1 vector of reaction indices
+0581 
+0582     % === 3) Group in one shot: unique subsystems + members per subsystem ===
+0583     [subSystems, ~, g] = unique(flatNames, 'stable');          % stable preserves first appearance
+0584     membersIdx = accumarray(g, flatIdx, [], @(v){v});          % cell: one vector of r-idx per group
+0585     nSubs = numel(subSystems);
+0586 
+0587     % === 4) Preallocate and build SBML groups (single simple loop) ===
+0588     modelSBML.groups_group(1:nSubs) = grpTemplate;
+0589     groupIDs = "group" + (1:nSubs);
+0590 
+0591     if isfield(grpTemplate,'groups_member')
+0592         memTemplate = grpTemplate.groups_member;
+0593     else
+0594         memTemplate = struct('groups_idRef','');
+0595     end
+0596 
+0597     for i = 1:nSubs
+0598         rIdx       = membersIdx{i};
+0599         groupRXNs  = model.rxns(rIdx);
+0600         cgroup     = grpTemplate;
+0601 
+0602         % Preallocate and fill members
+0603         nM = numel(groupRXNs);
+0604         if nM > 0
+0605             cgroup.groups_member(1:nM) = memTemplate;
+0606             [cgroup.groups_member.groups_idRef] = deal(groupRXNs{:});
+0607         else
+0608             cgroup.groups_member = repmat(memTemplate, 0, 1);
+0609         end
 0610 
-0611 modelSBML.fbc_objective.fbc_type='maximize';
-0612 modelSBML.fbc_objective.fbc_id='obj';
-0613 
-0614 ind=find(model.c);
-0615 
-0616 if isempty(ind)
-0617     modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0;
-0618 else
-0619     for i=1:length(ind)
-0620         %Copy the default values to the next index as long as it is not the
-0621         %last one
-0622         if i<numel(ind)
-0623             modelSBML.reaction(i+1)=modelSBML.reaction(i);
-0624         end
-0625         values=model.c(model.c~=0);
-0626         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id;
-0627         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i);
-0628         modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1;
-0629     end
-0630 end
-0631 
-0632 modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id;
-0633 
-0634 fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))];
-0635 if modelHasSubsystems
-0636     groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))];
-0637     modelSBML.namespaces=struct('prefix',{'','fbc','groups'},...
-0638     'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
-0639     fbcStr,groupStr});
-0640 else
-0641     modelSBML.namespaces=struct('prefix',{'','fbc'},...
-0642     'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
-0643     fbcStr});
-0644 end
-0645 
-0646 if sbmlPackageVersions(1) == 2
-0647     modelSBML.fbc_strict=1;
-0648     modelSBML.isSetfbc_strict = 1;
-0649 end
-0650 
-0651 modelSBML.rule=[];
-0652 modelSBML.constraint=[];
-0653 
-0654 [ravenDir,prevDir]=findRAVENroot();
-0655 fileName=checkFileExistence(fileName,1,true,false);
-0656 
-0657 OutputSBML_RAVEN(modelSBML,fileName,1,0,[1,0]);
+0611         cgroup.groups_id   = char(groupIDs(i));     % keep as char for SBML compatibility
+0612         cgroup.groups_name = subSystems{i};
+0613         modelSBML.groups_group(i) = cgroup;
+0614     end
+0615 end
+0616 
+0617 
+0618 %Prepare fbc_objective subfield
+0619 
+0620 modelSBML.fbc_objective.fbc_type='maximize';
+0621 modelSBML.fbc_objective.fbc_id='obj';
+0622 
+0623 ind=find(model.c);
+0624 
+0625 if isempty(ind)
+0626     modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0;
+0627 else
+0628     for i=1:length(ind)
+0629         %Copy the default values to the next index as long as it is not the
+0630         %last one
+0631         if i<numel(ind)
+0632             modelSBML.reaction(i+1)=modelSBML.reaction(i);
+0633         end
+0634         values=model.c(model.c~=0);
+0635         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id;
+0636         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i);
+0637         modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1;
+0638     end
+0639 end
+0640 
+0641 modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id;
+0642 
+0643 fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))];
+0644 if modelHasSubsystems
+0645     groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))];
+0646     modelSBML.namespaces=struct('prefix',{'','fbc','groups'},...
+0647         'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
+0648         fbcStr,groupStr});
+0649 else
+0650     modelSBML.namespaces=struct('prefix',{'','fbc'},...
+0651         'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
+0652         fbcStr});
+0653 end
+0654 
+0655 if sbmlPackageVersions(1) == 2
+0656     modelSBML.fbc_strict=1;
+0657     modelSBML.isSetfbc_strict = 1;
 0658 end
 0659 
-0660 
-0661 function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions)
-0662 %Returns the blank SBML model structure by using appropriate libSBML
-0663 %functions. This creates structure by considering three levels
-0664 
-0665 sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0666 sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0667 
-0668 for i=1:numel(sbmlFieldNames)
-0669     modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
-0670     sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0671     sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0672     if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment')
-0673         for j=1:numel(sbmlSubfieldNames)
-0674             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j};
-0675             sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0676             sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0677             if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw')
-0678                 for k=1:numel(sbmlSubsubfieldNames)
-0679                     %'compartment' and 'species' fields are not supposed to
-0680                     %have their standalone structures if they are subfields
-0681                     %or subsubfields
-0682                     if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species')
-0683                         modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k};
-0684                     end
-0685                     %If it is fbc_association in the third level, we need
-0686                     %to establish the fourth level, since libSBML requires
-0687                     %it
-0688                     if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association')
-0689                         fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0690                         fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0691                         for l=1:numel(fbc_associationFieldNames)
-0692                             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l};
-0693                         end
-0694                     end
-0695                 end
-0696             end
-0697         end
-0698     end
-0699     if ~isstruct(modelSBML.(sbmlFieldNames{1,i}))
-0700         modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
-0701     end
-0702 end
-0703 
-0704 modelSBML.unitDefinition.id='mmol_per_gDW_per_hr';
-0705 
-0706 unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0707 unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
-0708 
-0709 kinds={'mole','gram','second'};
-0710 exponents=[1 -1 -1];
-0711 scales=[-3 0 0];
-0712 multipliers=[1 1 1*60*60];
-0713 
-0714 for i=1:numel(unitFieldNames)
-0715     modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i};
-0716     for j=1:3
-0717         modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i};
-0718         if strcmp(unitFieldNames{1,i},'kind')
-0719             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j};
-0720         elseif strcmp(unitFieldNames{1,i},'exponent')
-0721             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j);
-0722         elseif strcmp(unitFieldNames{1,i},'scale')
-0723             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j);
-0724         elseif strcmp(unitFieldNames{1,i},'multiplier')
-0725             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j);
-0726         end
-0727     end
-0728 end
-0729 end
-0730 
-0731 function miriamString=getMiriam(miriamStruct)
-0732 %Returns a string with list elements for a miriam structure ('<rdf:li
-0733 %rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This
-0734 %is just to speed up things since this is done many times during the
-0735 %exporting
-0736 
-0737 miriamString='';
-0738 if isfield(miriamStruct,'name')
-0739     for i=1:numel(miriamStruct.name)
-0740         miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>'];
-0741     end
-0742 end
-0743 end
-0744 
-0745 function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i)
-0746 %This function provides reactants and products for particular reaction. The
-0747 %function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679
-0748 
-0749 met_idx = find(model.S(:, i));
-0750 tmp_Rxn.product=[];
-0751 tmp_Rxn.reactant=[];
-0752 for j_met=1:size(met_idx,1)
-0753     tmp_idx = met_idx(j_met,1);
-0754     sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id;
-0755     met_stoich = model.S(tmp_idx, i);
-0756     sbml_tmp_species_ref.stoichiometry = abs(met_stoich);
-0757     sbml_tmp_species_ref.isSetStoichiometry=1;
-0758     sbml_tmp_species_ref.constant=1;
-0759     if (met_stoich > 0)
-0760         tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ];
-0761     else
-0762         tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref];
-0763     end
-0764 end
-0765 end
-0766 
-0767 function vecT = columnVector(vec)
-0768 % Code below taken from COBRA Toolbox under GNU General Public License v3.0
-0769 % license file in readme/GPL.MD.
-0770 %
-0771 % Converts a vector to a column vector
-0772 %
-0773 % USAGE:
-0774 %
-0775 %   vecT = columnVector(vec)
-0776 %
-0777 % INPUT:
-0778 %   vec:     a vector
+0660 modelSBML.rule=[];
+0661 modelSBML.constraint=[];
+0662 
+0663 [ravenDir,prevDir]=findRAVENroot();
+0664 fileName=checkFileExistence(fileName,1,true,false);
+0665 
+0666 OutputSBML_RAVEN(modelSBML,fileName,1,0,[1,0]);
+0667 end
+0668 
+0669 
+0670 function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions)
+0671 %Returns the blank SBML model structure by using appropriate libSBML
+0672 %functions. This creates structure by considering three levels
+0673 
+0674 sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0675 sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0676 
+0677 for i=1:numel(sbmlFieldNames)
+0678     modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
+0679     sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0680     sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0681     if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment')
+0682         for j=1:numel(sbmlSubfieldNames)
+0683             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j};
+0684             sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0685             sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0686             if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw')
+0687                 for k=1:numel(sbmlSubsubfieldNames)
+0688                     %'compartment' and 'species' fields are not supposed to
+0689                     %have their standalone structures if they are subfields
+0690                     %or subsubfields
+0691                     if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species')
+0692                         modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k};
+0693                     end
+0694                     %If it is fbc_association in the third level, we need
+0695                     %to establish the fourth level, since libSBML requires
+0696                     %it
+0697                     if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association')
+0698                         fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0699                         fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0700                         for l=1:numel(fbc_associationFieldNames)
+0701                             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l};
+0702                         end
+0703                     end
+0704                 end
+0705             end
+0706         end
+0707     end
+0708     if ~isstruct(modelSBML.(sbmlFieldNames{1,i}))
+0709         modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
+0710     end
+0711 end
+0712 
+0713 modelSBML.unitDefinition.id='mmol_per_gDW_per_hr';
+0714 
+0715 unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0716 unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
+0717 
+0718 kinds={'mole','gram','second'};
+0719 exponents=[1 -1 -1];
+0720 scales=[-3 0 0];
+0721 multipliers=[1 1 1*60*60];
+0722 
+0723 for i=1:numel(unitFieldNames)
+0724     modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i};
+0725     for j=1:3
+0726         modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i};
+0727         if strcmp(unitFieldNames{1,i},'kind')
+0728             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j};
+0729         elseif strcmp(unitFieldNames{1,i},'exponent')
+0730             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j);
+0731         elseif strcmp(unitFieldNames{1,i},'scale')
+0732             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j);
+0733         elseif strcmp(unitFieldNames{1,i},'multiplier')
+0734             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j);
+0735         end
+0736     end
+0737 end
+0738 end
+0739 
+0740 function miriamString=getMiriam(miriamStruct)
+0741 %Returns a string with list elements for a miriam structure ('<rdf:li
+0742 %rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This
+0743 %is just to speed up things since this is done many times during the
+0744 %exporting
+0745 
+0746 miriamString='';
+0747 if isfield(miriamStruct,'name')
+0748     for i=1:numel(miriamStruct.name)
+0749         miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>'];
+0750     end
+0751 end
+0752 end
+0753 
+0754 function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i)
+0755 %This function provides reactants and products for particular reaction. The
+0756 %function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679
+0757 
+0758 met_idx = find(model.S(:, i));
+0759 tmp_Rxn.product=[];
+0760 tmp_Rxn.reactant=[];
+0761 for j_met=1:size(met_idx,1)
+0762     tmp_idx = met_idx(j_met,1);
+0763     sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id;
+0764     met_stoich = model.S(tmp_idx, i);
+0765     sbml_tmp_species_ref.stoichiometry = abs(met_stoich);
+0766     sbml_tmp_species_ref.isSetStoichiometry=1;
+0767     sbml_tmp_species_ref.constant=1;
+0768     if (met_stoich > 0)
+0769         tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ];
+0770     else
+0771         tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref];
+0772     end
+0773 end
+0774 end
+0775 
+0776 function vecT = columnVector(vec)
+0777 % Code below taken from COBRA Toolbox under GNU General Public License v3.0
+0778 % license file in readme/GPL.MD.
 0779 %
-0780 % OUTPUT:
-0781 %   vecT:    a column vector
-0782 
-0783 [n, m] = size(vec);
-0784 
-0785 if n < m
-0786     vecT = vec';
-0787 else
-0788     vecT = vec;
-0789 end
-0790 end
+0780 % Converts a vector to a column vector
+0781 %
+0782 % USAGE:
+0783 %
+0784 %   vecT = columnVector(vec)
+0785 %
+0786 % INPUT:
+0787 %   vec:     a vector
+0788 %
+0789 % OUTPUT:
+0790 %   vecT:    a column vector
+0791 
+0792 [n, m] = size(vec);
+0793 
+0794 if n < m
+0795     vecT = vec';
+0796 else
+0797     vecT = vec;
+0798 end
+0799 end
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file
diff --git a/io/exportModel.m b/io/exportModel.m
index 8beef3b2..ca969c17 100755
--- a/io/exportModel.m
+++ b/io/exportModel.m
@@ -245,7 +245,7 @@ function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
     if i'];
     end
-    
+
     % Export SBO terms from rxnMiriams
     if ~isempty(model.rxnMiriams{i})
         [~,sbo_ind] = ismember('sbo',model.rxnMiriams{i}.name);
@@ -491,7 +491,7 @@ function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
             model.rxnMiriams{i}.value(sbo_ind) = [];
         end
     end
-    
+
     %Export annotation information from rxnMiriams
     if (~isempty(model.rxnMiriams{i}) && isfield(modelSBML.reaction(i),'annotation')) || ~isempty(model.eccodes{i})
         modelSBML.reaction(i).annotation=[''];
@@ -504,14 +504,14 @@ function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
         end
         modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation getMiriam(model.rxnMiriams{i}) ''];
     end
-    
+
     if isfield(modelSBML.reaction, 'name')
         modelSBML.reaction(i).name=model.rxnNames{i};
     end
     if isfield(modelSBML.reaction, 'id')
         modelSBML.reaction(i).id=model.rxns{i};
     end
-    
+
     %Add the information about reactants and products
     involvedMets=addReactantsProducts(model,modelSBML,i);
     for j=1:numel(involvedMets.reactant)
@@ -553,59 +553,68 @@ function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
     modelSBML.reaction(i).fbc_upperFluxBound=totalNames{length(model.lb)+i};
 end
 
-%Prepare subSystems Code taken from COBRA functions getModelSubSystems,
-%writeSBML, findRxnsFromSubSystem under GNU General Public License v3.0,
-%license file in readme/GPL.MD. Code modified for RAVEN
+%Prepare subSystems
 if modelHasSubsystems
     modelSBML.groups_group.groups_kind = 'partonomy';
-    modelSBML.groups_group.sboTerm = 633;
-    tmpStruct=modelSBML.groups_group;
+    modelSBML.groups_group.sboTerm     = 633;
+    grpTemplate = modelSBML.groups_group;
 
-    rxns=model.rxns;
-    if ~any(cellfun(@iscell,model.subSystems))
-        if ~any(~cellfun(@isempty,model.subSystems))
-            subSystems = {};
-        else
-            subSystems = setdiff(model.subSystems,'');
-        end
+    % === 1) Normalize: make every entry a cell array of chars (vectorized) ===
+    isChar = cellfun(@ischar, model.subSystems);
+    model.subSystems(isChar)  = cellfun(@(s){s}, model.subSystems(isChar),  'UniformOutput', false);
+    model.subSystems(cellfun(@isempty, model.subSystems)) = {{}};
+
+    % If some entries contain string scalars by mistake, coerce them:
+    model.subSystems = cellfun(@(c) cellfun(@char, c, 'UniformOutput', false), model.subSystems, 'UniformOutput', false);
+
+    % === 2) Flatten once: names and their reaction indices (vectorized) ===
+    flatNames = [model.subSystems{:}];                    % 1×M cellstr of all subsystem labels
+    if isempty(flatNames)
+        % Nothing to do: no subsystems present
+        return
+    end
+
+    counts   = cellfun(@numel, model.subSystems);         % reactions -> how many subsystems
+    % For each reaction r, repeat r exactly counts(r) times
+    flatIdx  = arrayfun(@(r,c) repmat(r, c, 1), (1:numel(model.subSystems)).', counts, 'UniformOutput', false);
+    flatIdx  = vertcat(flatIdx{:});           % M×1 vector of reaction indices
+
+    % === 3) Group in one shot: unique subsystems + members per subsystem ===
+    [subSystems, ~, g] = unique(flatNames, 'stable');          % stable preserves first appearance
+    membersIdx = accumarray(g, flatIdx, [], @(v){v});          % cell: one vector of r-idx per group
+    nSubs = numel(subSystems);
+
+    % === 4) Preallocate and build SBML groups (single simple loop) ===
+    modelSBML.groups_group(1:nSubs) = grpTemplate;
+    groupIDs = "group" + (1:nSubs);
+
+    if isfield(grpTemplate,'groups_member')
+        memTemplate = grpTemplate.groups_member;
     else
-        orderedSubs = cellfun(@(x) columnVector(x),model.subSystems,'UniformOUtput',false);
-        subSystems = setdiff(vertcat(orderedSubs{:}),'');
-    end
-    if isempty(subSystems)
-        subSystems = {};
-    end
-    if ~isempty(subSystems)
-        %Build the groups for the group package
-        groupIDs = strcat('group',cellfun(@num2str, num2cell(1:length(subSystems)),'UniformOutput',false));
-        for i = 1:length(subSystems)
-            cgroup = tmpStruct;
-            if ~any(cellfun(@iscell,model.subSystems))
-                present = ismember(model.subSystems,subSystems{i});
-            else
-                present = cellfun(@(x) any(ismember(x,subSystems{i})),model.subSystems);
-            end
-            groupMembers = rxns(present);
-            for j = 1:numel(groupMembers)
-                cMember = tmpStruct.groups_member;
-                cMember.groups_idRef = groupMembers{j};
-                if j == 1
-                    cgroup.groups_member = cMember;
-                else
-                    cgroup.groups_member(j) = cMember;
-                end
-            end
-            cgroup.groups_id = groupIDs{i};
-            cgroup.groups_name = subSystems{i};
-            if i == 1
-                modelSBML.groups_group = cgroup;
-            else
-                modelSBML.groups_group(i) = cgroup;
-            end
+        memTemplate = struct('groups_idRef','');
+    end
+
+    for i = 1:nSubs
+        rIdx       = membersIdx{i};
+        groupRXNs  = model.rxns(rIdx);
+        cgroup     = grpTemplate;
+
+        % Preallocate and fill members
+        nM = numel(groupRXNs);
+        if nM > 0
+            cgroup.groups_member(1:nM) = memTemplate;
+            [cgroup.groups_member.groups_idRef] = deal(groupRXNs{:});
+        else
+            cgroup.groups_member = repmat(memTemplate, 0, 1);
         end
+
+        cgroup.groups_id   = char(groupIDs(i));     % keep as char for SBML compatibility
+        cgroup.groups_name = subSystems{i};
+        modelSBML.groups_group(i) = cgroup;
     end
 end
 
+
 %Prepare fbc_objective subfield
 
 modelSBML.fbc_objective.fbc_type='maximize';
@@ -635,12 +644,12 @@ function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
 if modelHasSubsystems
     groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))];
     modelSBML.namespaces=struct('prefix',{'','fbc','groups'},...
-    'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
-    fbcStr,groupStr});
+        'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
+        fbcStr,groupStr});
 else
     modelSBML.namespaces=struct('prefix',{'','fbc'},...
-    'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
-    fbcStr});
+        'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
+        fbcStr});
 end
 
 if sbmlPackageVersions(1) == 2

From 9fd7e4b9d120a34ec75e17b29b79cda469803980 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven 
Date: Thu, 15 Jan 2026 22:50:08 +0100
Subject: [PATCH 9/9] feat: runRAVENtests for quick local tests

---
 doc/index.html               | 68 ++++++++++++++++++------------------
 doc/utils/index.html         |  2 +-
 doc/utils/runRAVENtests.html | 60 +++++++++++++++++++++++++++++++
 utils/runRAVENtests.m        | 11 ++++++
 4 files changed, 106 insertions(+), 35 deletions(-)
 create mode 100644 doc/utils/runRAVENtests.html
 create mode 100644 utils/runRAVENtests.m

diff --git a/doc/index.html b/doc/index.html
index faf716ce..025e7fcb 100644
--- a/doc/index.html
+++ b/doc/index.html
@@ -38,40 +38,40 @@ 
    Matlab Files found in these Directories
    
 				buildEquation		fitTasks		importExcelModel		runINIT		
 				canConsume		followChanged		importExportTests		runPhenotypePhasePlane		
 				canProduce		followFluxes		importModel		runProductionEnvelope		
-				cdhitTests		ftINIT		linkMetaCycKEGGRxns		runRobustnessAnalysis		
-				changeGeneAssoc		ftINITFillGaps		loadSheet		runSimpleOptKnock		
-				changeGrRules		ftINITFillGapsForAllTasks		loadWorkbook		scoreComplexModel		
-				changeRxns		ftINITFillGapsMILP		mafftTests		scoreModel		
-				checkFileExistence		ftINITInternalAlg		makeFakeBlastStructure		setColorToMapRxns		
-				checkFunctionUniqueness		gapReport		makeSomething		setExchangeBounds		
-				checkInstallation		generateNewIds		mapCompartments		setOmicDataToRxns		
-				checkModelStruct		getAllRxnsFromGenes		mapPathwayRxnNames		setParam		
-				checkProduction		getAllSubGraphs		markPathwayWithExpression		setRavenSolver		
-				checkRxn		getAllowedBounds		markPathwayWithFluxes		setTitle		
-				checkSolution		getBlast		mergeCompartments		simplifyModel		
-				checkTasks		getBlastFromExcel		mergeLinear		solveLP		
-				checkTasksTests		getColorCodes		mergeModels		solveQP		
-				cleanSheet		getDiamond		miriamTests		solverTests		
-				closeModel		getElementalBalance		modelAbilitiesTests		sortIdentifiers		
-				colorPathway		getEnzymesFromMetaCyc		modelConversionTests		sortModel		
-				colorSubsystem		getEssentialRxns		modelCurationTests		standardizeGrRules		
-				combineMetaCycKEGGModels		getExchangeRxns		modelSortingTests		standardizeModelFieldOrder		
-				compareMultipleModels		getExprForRxnScore		optimizeProb		startup		
-				compareRxnsGenesMetsComps		getExpressionStructure		parallelPoolRAVEN		tinitTests		
-				constructEquations		getFluxZ		parseFormulas		trimPathway		
-				constructMultiFasta		getFullPath		parseHPA		tutorial1		
-				constructPathwayFromCelldesigner		getGenesFromGrRules		parseHPArna		tutorial2		
-				constructS		getGenesFromKEGG		parseRxnEqu		tutorial2_solutions		
-				consumeSomething		getINITModel		parseScores		tutorial3		
-				contractModel		getINITSteps		parseTaskList		tutorial3_solutions		
-				convertCharArray		getIndexes		permuteModel		tutorial4		
-				convertToIrrev		getKEGGModelForOrganism		plotAdditionalInfo		tutorial4_solutions		
-				copyToComps		getMD5Hash		plotLabels		tutorial5		
-				deleteUnusedGenes		getMetaCycModelForOrganism		predictLocalization		tutorial6		
-				diamondTests		getMetsFromKEGG		prepINITModel		updateDocumentation		
-				dispEM		getMetsFromMetaCyc		printFluxes		writeSheet		
-				drawMap		getMetsInComp		printModel		writeYAMLmodel		
-				drawPathway		getMinNrFluxes		printModelStats				
+				cdhitTests		ftINIT		linkMetaCycKEGGRxns		runRAVENtests		
+				changeGeneAssoc		ftINITFillGaps		loadSheet		runRobustnessAnalysis		
+				changeGrRules		ftINITFillGapsForAllTasks		loadWorkbook		runSimpleOptKnock		
+				changeRxns		ftINITFillGapsMILP		mafftTests		scoreComplexModel		
+				checkFileExistence		ftINITInternalAlg		makeFakeBlastStructure		scoreModel		
+				checkFunctionUniqueness		gapReport		makeSomething		setColorToMapRxns		
+				checkInstallation		generateNewIds		mapCompartments		setExchangeBounds		
+				checkModelStruct		getAllRxnsFromGenes		mapPathwayRxnNames		setOmicDataToRxns		
+				checkProduction		getAllSubGraphs		markPathwayWithExpression		setParam		
+				checkRxn		getAllowedBounds		markPathwayWithFluxes		setRavenSolver		
+				checkSolution		getBlast		mergeCompartments		setTitle		
+				checkTasks		getBlastFromExcel		mergeLinear		simplifyModel		
+				checkTasksTests		getColorCodes		mergeModels		solveLP		
+				cleanSheet		getDiamond		miriamTests		solveQP		
+				closeModel		getElementalBalance		modelAbilitiesTests		solverTests		
+				colorPathway		getEnzymesFromMetaCyc		modelConversionTests		sortIdentifiers		
+				colorSubsystem		getEssentialRxns		modelCurationTests		sortModel		
+				combineMetaCycKEGGModels		getExchangeRxns		modelSortingTests		standardizeGrRules		
+				compareMultipleModels		getExprForRxnScore		optimizeProb		standardizeModelFieldOrder		
+				compareRxnsGenesMetsComps		getExpressionStructure		parallelPoolRAVEN		startup		
+				constructEquations		getFluxZ		parseFormulas		tinitTests		
+				constructMultiFasta		getFullPath		parseHPA		trimPathway		
+				constructPathwayFromCelldesigner		getGenesFromGrRules		parseHPArna		tutorial1		
+				constructS		getGenesFromKEGG		parseRxnEqu		tutorial2		
+				consumeSomething		getINITModel		parseScores		tutorial2_solutions		
+				contractModel		getINITSteps		parseTaskList		tutorial3		
+				convertCharArray		getIndexes		permuteModel		tutorial3_solutions		
+				convertToIrrev		getKEGGModelForOrganism		plotAdditionalInfo		tutorial4		
+				copyToComps		getMD5Hash		plotLabels		tutorial4_solutions		
+				deleteUnusedGenes		getMetaCycModelForOrganism		predictLocalization		tutorial5		
+				diamondTests		getMetsFromKEGG		prepINITModel		tutorial6		
+				dispEM		getMetsFromMetaCyc		printFluxes		updateDocumentation		
+				drawMap		getMetsInComp		printModel		writeSheet		
+				drawPathway		getMinNrFluxes		printModelStats		writeYAMLmodel		
 				editMiriam		getModelFromHomology		printOrange			
 
 
diff --git a/doc/utils/index.html b/doc/utils/index.html
index f71a3bf0..9a1637b6 100644
--- a/doc/utils/index.html
+++ b/doc/utils/index.html
@@ -19,7 +19,7 @@ 
    Index for utils
    
 
 
    Matlab files in this directory:
    
 
-
     emptyOrLogicalScalar 
     emptyOrTextOrCellOfTextValidate [] OR text scalar (char row or string scalar) OR cell array of such text. 
     emptyOrTextScalar 
    
+ emptyOrLogicalScalar  emptyOrTextOrCellOfTextValidate [] OR text scalar (char row or string scalar) OR cell array of such text.  emptyOrTextScalar  runRAVENtestsrunRAVENtests 
 
 
 
diff --git a/doc/utils/runRAVENtests.html b/doc/utils/runRAVENtests.html
new file mode 100644
index 00000000..da9c79b7
--- /dev/null
+++ b/doc/utils/runRAVENtests.html
@@ -0,0 +1,60 @@
+
+
+
+  Description of runRAVENtests
+  
+  
+  
+  
+  
+  
+
+
+
+
    Home >  utils > runRAVENtests.m
    
+
+
+
+
    runRAVENtests
+
    
+
+
    PURPOSE ^
    
+
    runRAVENtests
    
+
+
    SYNOPSIS ^
    
+
    function testResults = runRAVENtests 
    
+
+
    DESCRIPTION ^
    
+
     runRAVENtests
+  Runs all unit tests found at RAVEN/test/unit_tests, and prints output in
+  command window.
    
+
+
+
    CROSS-REFERENCE INFORMATION ^
    
+This function calls:
+
      
+
    
+This function is called by:
+
      
+
    
+
+
+
+
+
    SOURCE CODE ^
    
+
    0001 function testResults = runRAVENtests
+0002 % runRAVENtests
+0003 %  Runs all unit tests found at RAVEN/test/unit_tests, and prints output in
+0004 %  command window.
+0005 
+0006 ravenPath = findRAVENroot;
+0007 curwd = pwd;
+0008 cd(fullfile(ravenPath,'testing','unit_tests'));
+0009 testResults = runtests(struct2table(dir('*.m')).name);
+0010 
+0011 cd(curwd);
    
+
    Generated by m2html © 2005
    
+
+
\ No newline at end of file
diff --git a/utils/runRAVENtests.m b/utils/runRAVENtests.m
new file mode 100644
index 00000000..cd72cc6a
--- /dev/null
+++ b/utils/runRAVENtests.m
@@ -0,0 +1,11 @@
+function testResults = runRAVENtests
+% runRAVENtests
+%  Runs all unit tests found at RAVEN/test/unit_tests, and prints output in
+%  command window.
+
+ravenPath = findRAVENroot;
+curwd = pwd;
+cd(fullfile(ravenPath,'testing','unit_tests'));
+testResults = runtests(struct2table(dir('*.m')).name);
+
+cd(curwd);
\ No newline at end of file