From 62e7546664b5d24ede8cd801612cd1c8f7b5434a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 28 Oct 2024 21:27:24 +0100 Subject: [PATCH 1/9] fix: addMets if metNames and not mets are given --- core/addMets.m | 1 + doc/core/addMets.html | 597 +++++++++++++++++++++--------------------- 2 files changed, 300 insertions(+), 298 deletions(-) diff --git a/core/addMets.m b/core/addMets.m index 57bceef8..3b4668a6 100755 --- a/core/addMets.m +++ b/core/addMets.m @@ -70,6 +70,7 @@ %Check some stuff regarding the required fields if ~isfield(metsToAdd,'mets') + metsToAdd.metNames=convertCharArray(metsToAdd.metNames); metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames)); else metsToAdd.mets=convertCharArray(metsToAdd.mets); diff --git a/doc/core/addMets.html b/doc/core/addMets.html index ff66f36e..da17e221 100644 --- a/doc/core/addMets.html +++ b/doc/core/addMets.html @@ -166,306 +166,307 @@

SOURCE CODE ^%Check some stuff regarding the required fields 0072 if ~isfield(metsToAdd,'mets') -0073 metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames)); -0074 else -0075 metsToAdd.mets=convertCharArray(metsToAdd.mets); -0076 end -0077 if ~isfield(metsToAdd,'metNames') -0078 metsToAdd.metNames=metsToAdd.mets; -0079 else -0080 metsToAdd.metNames=convertCharArray(metsToAdd.metNames); -0081 end -0082 if ~isfield(metsToAdd,'compartments') -0083 EM='compartments is a required field in metsToAdd'; -0084 dispEM(EM); -0085 else -0086 metsToAdd.compartments=convertCharArray(metsToAdd.compartments); -0087 %If only one compartment is given, assume it is for all metabolites -0088 if numel(metsToAdd.compartments)==1 && numel(metsToAdd.mets)>1 -0089 temp=cell(numel(metsToAdd.mets),1); -0090 temp(:)=metsToAdd.compartments; -0091 metsToAdd.compartments=temp; -0092 end -0093 end -0094 -0095 %Number of metabolites -0096 nMets=numel(metsToAdd.mets); -0097 nOldMets=numel(model.mets); -0098 filler=cell(nMets,1); -0099 filler(:)={''}; -0100 largeFiller=cell(nOldMets,1); -0101 largeFiller(:)={''}; -0102 -0103 %Check that no metabolite ids are already present in the model -0104 I=ismember(metsToAdd.mets,model.mets); -0105 if any(I) -0106 error('One or more elements in metsToAdd.mets are already present in model.mets: %s',... -0107 strjoin(metsToAdd.mets(I),', ')); -0108 else -0109 newModel.mets=[newModel.mets;metsToAdd.mets(:)]; -0110 end -0111 -0112 %Check that all the compartments could be found -0113 [I, compMap]=ismember(metsToAdd.compartments,model.comps); -0114 if ~all(I) -0115 EM='metsToAdd.compartments must match model.comps'; -0116 dispEM(EM); -0117 end -0118 -0119 %Check that the metabolite names aren't present in the same compartment. -0120 %Not the neatest way maybe.. -0121 t1=strcat(metsToAdd.metNames(:),'***',metsToAdd.compartments(:)); -0122 t2=strcat(model.metNames,'***',model.comps(model.metComps)); -0123 if any(ismember(t1,t2)) -0124 EM='One or more elements in metsToAdd.metNames already exist in the same compartments as the one it is being added to'; -0125 dispEM(EM); -0126 end -0127 -0128 %Some more checks and if they pass then add each field to the structure -0129 if numel(metsToAdd.metNames)~=nMets -0130 EM='metsToAdd.metNames must have the same number of elements as metsToAdd.mets'; -0131 dispEM(EM); -0132 else -0133 newModel.metNames=[newModel.metNames;metsToAdd.metNames(:)]; -0134 end -0135 -0136 if numel(compMap)~=nMets -0137 EM='metsToAdd.compartments must have the same number of elements as metsToAdd.mets'; -0138 dispEM(EM); -0139 else -0140 newModel.metComps=[newModel.metComps;compMap(:)]; -0141 end -0142 -0143 if isfield(metsToAdd,'b') -0144 if size(metsToAdd.b,1)~=nMets -0145 EM='metsToAdd.b must have the same number of elements as metsToAdd.mets'; -0146 dispEM(EM); -0147 else -0148 %Add empty field if it doesn't exist -0149 if ~isfield(newModel,'b') -0150 newModel.b=zeros(nOldMets,1); -0151 end -0152 -0153 %If the original is only one vector -0154 if size(metsToAdd.b,2)>size(newModel.b,2) -0155 newModel.b=[newModel.b newModel.b]; -0156 end -0157 %Add the new ones -0158 newModel.b=[newModel.b;metsToAdd.b]; -0159 end -0160 else -0161 if isfield(newModel,'b') -0162 %Add the default -0163 newModel.b=[newModel.b;zeros(nMets,size(newModel.b,2))]; -0164 end -0165 end -0166 -0167 if isfield(metsToAdd,'unconstrained') -0168 if numel(metsToAdd.unconstrained)~=nMets -0169 EM='metsToAdd.unconstrained must have the same number of elements as metsToAdd.mets'; -0170 dispEM(EM); -0171 else -0172 %Add empty field if it doesn't exist -0173 if ~isfield(newModel,'unconstrained') -0174 newModel.unconstrained=zeros(nOldMets,1); -0175 end -0176 -0177 %Add the new ones -0178 newModel.unconstrained=[newModel.unconstrained;metsToAdd.unconstrained(:)]; -0179 end -0180 else -0181 if isfield(newModel,'unconstrained') -0182 %Add the default -0183 newModel.unconstrained=[newModel.unconstrained;zeros(nMets,1)]; -0184 end -0185 end -0186 -0187 if isfield(metsToAdd,'inchis') -0188 metsToAdd.inchis=convertCharArray(metsToAdd.inchis); -0189 if numel(metsToAdd.inchis)~=nMets -0190 EM='metsToAdd.inchis must have the same number of elements as metsToAdd.mets'; -0191 dispEM(EM); -0192 end -0193 %Add empty field if it doesn't exist -0194 if ~isfield(newModel,'inchis') -0195 newModel.inchis=largeFiller; -0196 end -0197 newModel.inchis=[newModel.inchis;metsToAdd.inchis(:)]; -0198 else -0199 %Add empty strings if structure is in model -0200 if isfield(newModel,'inchis') -0201 newModel.inchis=[newModel.inchis;filler]; -0202 end -0203 end -0204 +0073 metsToAdd.metNames=convertCharArray(metsToAdd.metNames); +0074 metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames)); +0075 else +0076 metsToAdd.mets=convertCharArray(metsToAdd.mets); +0077 end +0078 if ~isfield(metsToAdd,'metNames') +0079 metsToAdd.metNames=metsToAdd.mets; +0080 else +0081 metsToAdd.metNames=convertCharArray(metsToAdd.metNames); +0082 end +0083 if ~isfield(metsToAdd,'compartments') +0084 EM='compartments is a required field in metsToAdd'; +0085 dispEM(EM); +0086 else +0087 metsToAdd.compartments=convertCharArray(metsToAdd.compartments); +0088 %If only one compartment is given, assume it is for all metabolites +0089 if numel(metsToAdd.compartments)==1 && numel(metsToAdd.mets)>1 +0090 temp=cell(numel(metsToAdd.mets),1); +0091 temp(:)=metsToAdd.compartments; +0092 metsToAdd.compartments=temp; +0093 end +0094 end +0095 +0096 %Number of metabolites +0097 nMets=numel(metsToAdd.mets); +0098 nOldMets=numel(model.mets); +0099 filler=cell(nMets,1); +0100 filler(:)={''}; +0101 largeFiller=cell(nOldMets,1); +0102 largeFiller(:)={''}; +0103 +0104 %Check that no metabolite ids are already present in the model +0105 I=ismember(metsToAdd.mets,model.mets); +0106 if any(I) +0107 error('One or more elements in metsToAdd.mets are already present in model.mets: %s',... +0108 strjoin(metsToAdd.mets(I),', ')); +0109 else +0110 newModel.mets=[newModel.mets;metsToAdd.mets(:)]; +0111 end +0112 +0113 %Check that all the compartments could be found +0114 [I, compMap]=ismember(metsToAdd.compartments,model.comps); +0115 if ~all(I) +0116 EM='metsToAdd.compartments must match model.comps'; +0117 dispEM(EM); +0118 end +0119 +0120 %Check that the metabolite names aren't present in the same compartment. +0121 %Not the neatest way maybe.. +0122 t1=strcat(metsToAdd.metNames(:),'***',metsToAdd.compartments(:)); +0123 t2=strcat(model.metNames,'***',model.comps(model.metComps)); +0124 if any(ismember(t1,t2)) +0125 EM='One or more elements in metsToAdd.metNames already exist in the same compartments as the one it is being added to'; +0126 dispEM(EM); +0127 end +0128 +0129 %Some more checks and if they pass then add each field to the structure +0130 if numel(metsToAdd.metNames)~=nMets +0131 EM='metsToAdd.metNames must have the same number of elements as metsToAdd.mets'; +0132 dispEM(EM); +0133 else +0134 newModel.metNames=[newModel.metNames;metsToAdd.metNames(:)]; +0135 end +0136 +0137 if numel(compMap)~=nMets +0138 EM='metsToAdd.compartments must have the same number of elements as metsToAdd.mets'; +0139 dispEM(EM); +0140 else +0141 newModel.metComps=[newModel.metComps;compMap(:)]; +0142 end +0143 +0144 if isfield(metsToAdd,'b') +0145 if size(metsToAdd.b,1)~=nMets +0146 EM='metsToAdd.b must have the same number of elements as metsToAdd.mets'; +0147 dispEM(EM); +0148 else +0149 %Add empty field if it doesn't exist +0150 if ~isfield(newModel,'b') +0151 newModel.b=zeros(nOldMets,1); +0152 end +0153 +0154 %If the original is only one vector +0155 if size(metsToAdd.b,2)>size(newModel.b,2) +0156 newModel.b=[newModel.b newModel.b]; +0157 end +0158 %Add the new ones +0159 newModel.b=[newModel.b;metsToAdd.b]; +0160 end +0161 else +0162 if isfield(newModel,'b') +0163 %Add the default +0164 newModel.b=[newModel.b;zeros(nMets,size(newModel.b,2))]; +0165 end +0166 end +0167 +0168 if isfield(metsToAdd,'unconstrained') +0169 if numel(metsToAdd.unconstrained)~=nMets +0170 EM='metsToAdd.unconstrained must have the same number of elements as metsToAdd.mets'; +0171 dispEM(EM); +0172 else +0173 %Add empty field if it doesn't exist +0174 if ~isfield(newModel,'unconstrained') +0175 newModel.unconstrained=zeros(nOldMets,1); +0176 end +0177 +0178 %Add the new ones +0179 newModel.unconstrained=[newModel.unconstrained;metsToAdd.unconstrained(:)]; +0180 end +0181 else +0182 if isfield(newModel,'unconstrained') +0183 %Add the default +0184 newModel.unconstrained=[newModel.unconstrained;zeros(nMets,1)]; +0185 end +0186 end +0187 +0188 if isfield(metsToAdd,'inchis') +0189 metsToAdd.inchis=convertCharArray(metsToAdd.inchis); +0190 if numel(metsToAdd.inchis)~=nMets +0191 EM='metsToAdd.inchis must have the same number of elements as metsToAdd.mets'; +0192 dispEM(EM); +0193 end +0194 %Add empty field if it doesn't exist +0195 if ~isfield(newModel,'inchis') +0196 newModel.inchis=largeFiller; +0197 end +0198 newModel.inchis=[newModel.inchis;metsToAdd.inchis(:)]; +0199 else +0200 %Add empty strings if structure is in model +0201 if isfield(newModel,'inchis') +0202 newModel.inchis=[newModel.inchis;filler]; +0203 end +0204 end 0205 -0206 if isfield(metsToAdd,'metSmiles') -0207 metsToAdd.metSmiles=convertCharArray(metsToAdd.metSmiles); -0208 if numel(metsToAdd.metSmiles)~=nMets -0209 EM='metsToAdd.metSmiles must have the same number of elements as metsToAdd.mets'; -0210 dispEM(EM); -0211 end -0212 %Add empty field if it doesn't exist -0213 if ~isfield(newModel,'metSmiles') -0214 newModel.metSmiles=largeFiller; -0215 end -0216 newModel.metSmiles=[newModel.metSmiles;metsToAdd.metSmiles(:)]; -0217 else -0218 %Add empty strings if structure is in model -0219 if isfield(newModel,'metSmiles') -0220 newModel.metSmiles=[newModel.metSmiles;filler]; -0221 end -0222 end -0223 -0224 if isfield(metsToAdd,'metFormulas') -0225 metsToAdd.metFormulas=convertCharArray(metsToAdd.metFormulas); -0226 if numel(metsToAdd.metFormulas)~=nMets -0227 EM='metsToAdd.metFormulas must have the same number of elements as metsToAdd.mets'; -0228 dispEM(EM); -0229 end -0230 %Add empty field if it doesn't exist -0231 if ~isfield(newModel,'metFormulas') -0232 newModel.metFormulas=largeFiller; -0233 end -0234 newModel.metFormulas=[newModel.metFormulas;metsToAdd.metFormulas(:)]; -0235 else -0236 %Add default -0237 if isfield(newModel,'metFormulas') -0238 newModel.metFormulas=[newModel.metFormulas;filler]; -0239 end -0240 end -0241 -0242 if isfield(metsToAdd,'metCharges') -0243 if numel(metsToAdd.metCharges)~=nMets -0244 EM='metsToAdd.metCharges must have the same number of elements as metsToAdd.mets'; -0245 dispEM(EM); -0246 end -0247 if ~isnumeric(metsToAdd.metCharges) -0248 EM='metsToAdd.metCharges must be of type "double"'; -0249 dispEM(EM); -0250 end -0251 if ~isfield(newModel,'metCharges') -0252 newModel.metCharges=NaN(numel(largeFiller),1); -0253 end -0254 newModel.metCharges=[newModel.metCharges;metsToAdd.metCharges(:)]; -0255 else -0256 %Add default -0257 if isfield(newModel,'metCharges') -0258 newModel.metCharges=[newModel.metCharges;NaN(numel(filler),1)]; -0259 end -0260 end -0261 -0262 if isfield(metsToAdd,'metDeltaG') -0263 if numel(metsToAdd.metDeltaG)~=nMets -0264 EM='metsToAdd.metDeltaG must have the same number of elements as metsToAdd.mets'; -0265 dispEM(EM); -0266 end -0267 if ~isnumeric(metsToAdd.metDeltaG) -0268 EM='metsToAdd.metDeltaG must be of type "double"'; -0269 dispEM(EM); -0270 end -0271 if ~isfield(newModel,'metDeltaG') -0272 newModel.metDeltaG=NaN(numel(largeFiller),1); -0273 end -0274 newModel.metDeltaG=[newModel.metDeltaG;metsToAdd.metDeltaG(:)]; -0275 else -0276 %Add default -0277 if isfield(newModel,'metDeltaG') -0278 newModel.metDeltaG=[newModel.metDeltaG;NaN(numel(filler),1)]; -0279 end -0280 end -0281 +0206 +0207 if isfield(metsToAdd,'metSmiles') +0208 metsToAdd.metSmiles=convertCharArray(metsToAdd.metSmiles); +0209 if numel(metsToAdd.metSmiles)~=nMets +0210 EM='metsToAdd.metSmiles must have the same number of elements as metsToAdd.mets'; +0211 dispEM(EM); +0212 end +0213 %Add empty field if it doesn't exist +0214 if ~isfield(newModel,'metSmiles') +0215 newModel.metSmiles=largeFiller; +0216 end +0217 newModel.metSmiles=[newModel.metSmiles;metsToAdd.metSmiles(:)]; +0218 else +0219 %Add empty strings if structure is in model +0220 if isfield(newModel,'metSmiles') +0221 newModel.metSmiles=[newModel.metSmiles;filler]; +0222 end +0223 end +0224 +0225 if isfield(metsToAdd,'metFormulas') +0226 metsToAdd.metFormulas=convertCharArray(metsToAdd.metFormulas); +0227 if numel(metsToAdd.metFormulas)~=nMets +0228 EM='metsToAdd.metFormulas must have the same number of elements as metsToAdd.mets'; +0229 dispEM(EM); +0230 end +0231 %Add empty field if it doesn't exist +0232 if ~isfield(newModel,'metFormulas') +0233 newModel.metFormulas=largeFiller; +0234 end +0235 newModel.metFormulas=[newModel.metFormulas;metsToAdd.metFormulas(:)]; +0236 else +0237 %Add default +0238 if isfield(newModel,'metFormulas') +0239 newModel.metFormulas=[newModel.metFormulas;filler]; +0240 end +0241 end +0242 +0243 if isfield(metsToAdd,'metCharges') +0244 if numel(metsToAdd.metCharges)~=nMets +0245 EM='metsToAdd.metCharges must have the same number of elements as metsToAdd.mets'; +0246 dispEM(EM); +0247 end +0248 if ~isnumeric(metsToAdd.metCharges) +0249 EM='metsToAdd.metCharges must be of type "double"'; +0250 dispEM(EM); +0251 end +0252 if ~isfield(newModel,'metCharges') +0253 newModel.metCharges=NaN(numel(largeFiller),1); +0254 end +0255 newModel.metCharges=[newModel.metCharges;metsToAdd.metCharges(:)]; +0256 else +0257 %Add default +0258 if isfield(newModel,'metCharges') +0259 newModel.metCharges=[newModel.metCharges;NaN(numel(filler),1)]; +0260 end +0261 end +0262 +0263 if isfield(metsToAdd,'metDeltaG') +0264 if numel(metsToAdd.metDeltaG)~=nMets +0265 EM='metsToAdd.metDeltaG must have the same number of elements as metsToAdd.mets'; +0266 dispEM(EM); +0267 end +0268 if ~isnumeric(metsToAdd.metDeltaG) +0269 EM='metsToAdd.metDeltaG must be of type "double"'; +0270 dispEM(EM); +0271 end +0272 if ~isfield(newModel,'metDeltaG') +0273 newModel.metDeltaG=NaN(numel(largeFiller),1); +0274 end +0275 newModel.metDeltaG=[newModel.metDeltaG;metsToAdd.metDeltaG(:)]; +0276 else +0277 %Add default +0278 if isfield(newModel,'metDeltaG') +0279 newModel.metDeltaG=[newModel.metDeltaG;NaN(numel(filler),1)]; +0280 end +0281 end 0282 -0283 if isfield(metsToAdd,'metNotes') -0284 metsToAdd.metNotes=convertCharArray(metsToAdd.metNotes); -0285 if numel(metsToAdd.metNotes)==1 && numel(metsToAdd.mets)>1 -0286 temp=cell(numel(metsToAdd.mets),1); -0287 temp(:)=metsToAdd.metNotes; -0288 metsToAdd.metNotes=temp; -0289 end -0290 if numel(metsToAdd.metNotes)~=nMets -0291 EM='metsToAdd.metNotes must have the same number of elements as metsToAdd.mets'; -0292 dispEM(EM); -0293 end -0294 %Add empty field if it doesn't exist -0295 if ~isfield(newModel,'metNotes') -0296 newModel.metNotes=largeFiller; -0297 end -0298 newModel.metNotes=[newModel.metNotes;metsToAdd.metNotes(:)]; -0299 else -0300 %Add empty strings if structure is in model -0301 if isfield(newModel,'metNotes') -0302 newModel.metNotes=[newModel.metNotes;filler]; -0303 end -0304 end -0305 -0306 %Don't check the type of metMiriams -0307 if isfield(metsToAdd,'metMiriams') -0308 if numel(metsToAdd.metMiriams)==1 && numel(metsToAdd.mets)>1 -0309 temp=cell(numel(metsToAdd.mets),1); -0310 temp(:)={metsToAdd.metMiriams}; -0311 metsToAdd.metMiriams=temp; -0312 end -0313 if numel(metsToAdd.metMiriams)~=nMets -0314 EM='metsToAdd.metMiriams must have the same number of elements as metsToAdd.mets'; -0315 dispEM(EM); -0316 end -0317 %Add empty field if it doesn't exist -0318 if ~isfield(newModel,'metMiriams') -0319 newModel.metMiriams=cell(nOldMets,1); -0320 end -0321 newModel.metMiriams=[newModel.metMiriams;metsToAdd.metMiriams(:)]; -0322 else -0323 if isfield(newModel,'metMiriams') -0324 newModel.metMiriams=[newModel.metMiriams;cell(nMets,1)]; -0325 end -0326 end -0327 -0328 if isfield(newModel,'metFrom') -0329 newModel.metFrom=[newModel.metFrom;filler]; -0330 end -0331 -0332 %Expand the S matrix -0333 newModel.S=[newModel.S;sparse(nMets,size(newModel.S,2))]; -0334 -0335 if copyInfo==true -0336 [I, J]=ismember(metsToAdd.metNames,model.metNames); -0337 J=J(I); -0338 %I is the indexes of the new metabolites for which a metabolite with -0339 %the same name existed -0340 I=find(I)+nOldMets; -0341 %Go through each of the added mets and copy annotation if it doesn't -0342 %exist -0343 for i=1:numel(I) -0344 if isfield(newModel,'inchis') -0345 if isempty(newModel.inchis{I(i)}) -0346 newModel.inchis(I(i))=newModel.inchis(J(i)); -0347 end -0348 end -0349 if isfield(newModel,'metSmiles') -0350 if isempty(newModel.metSmiles{I(i)}) -0351 newModel.metSmiles(I(i))=newModel.metSmiles(J(i)); -0352 end -0353 end -0354 if isfield(newModel,'metFormulas') -0355 if isempty(newModel.metFormulas{I(i)}) -0356 newModel.metFormulas(I(i))=newModel.metFormulas(J(i)); -0357 end -0358 end -0359 if isfield(newModel,'metMiriams') -0360 if isempty(newModel.metMiriams{I(i)}) -0361 newModel.metMiriams(I(i))=newModel.metMiriams(J(i)); -0362 end -0363 end -0364 if isfield(newModel,'metCharges') -0365 newModel.metCharges(I(i))=newModel.metCharges(J(i)); -0366 end -0367 if isfield(newModel,'metDeltaG') -0368 newModel.metDeltaG(I(i))=newModel.metDeltaG(J(i)); -0369 end -0370 end -0371 end -0372 end +0283 +0284 if isfield(metsToAdd,'metNotes') +0285 metsToAdd.metNotes=convertCharArray(metsToAdd.metNotes); +0286 if numel(metsToAdd.metNotes)==1 && numel(metsToAdd.mets)>1 +0287 temp=cell(numel(metsToAdd.mets),1); +0288 temp(:)=metsToAdd.metNotes; +0289 metsToAdd.metNotes=temp; +0290 end +0291 if numel(metsToAdd.metNotes)~=nMets +0292 EM='metsToAdd.metNotes must have the same number of elements as metsToAdd.mets'; +0293 dispEM(EM); +0294 end +0295 %Add empty field if it doesn't exist +0296 if ~isfield(newModel,'metNotes') +0297 newModel.metNotes=largeFiller; +0298 end +0299 newModel.metNotes=[newModel.metNotes;metsToAdd.metNotes(:)]; +0300 else +0301 %Add empty strings if structure is in model +0302 if isfield(newModel,'metNotes') +0303 newModel.metNotes=[newModel.metNotes;filler]; +0304 end +0305 end +0306 +0307 %Don't check the type of metMiriams +0308 if isfield(metsToAdd,'metMiriams') +0309 if numel(metsToAdd.metMiriams)==1 && numel(metsToAdd.mets)>1 +0310 temp=cell(numel(metsToAdd.mets),1); +0311 temp(:)={metsToAdd.metMiriams}; +0312 metsToAdd.metMiriams=temp; +0313 end +0314 if numel(metsToAdd.metMiriams)~=nMets +0315 EM='metsToAdd.metMiriams must have the same number of elements as metsToAdd.mets'; +0316 dispEM(EM); +0317 end +0318 %Add empty field if it doesn't exist +0319 if ~isfield(newModel,'metMiriams') +0320 newModel.metMiriams=cell(nOldMets,1); +0321 end +0322 newModel.metMiriams=[newModel.metMiriams;metsToAdd.metMiriams(:)]; +0323 else +0324 if isfield(newModel,'metMiriams') +0325 newModel.metMiriams=[newModel.metMiriams;cell(nMets,1)]; +0326 end +0327 end +0328 +0329 if isfield(newModel,'metFrom') +0330 newModel.metFrom=[newModel.metFrom;filler]; +0331 end +0332 +0333 %Expand the S matrix +0334 newModel.S=[newModel.S;sparse(nMets,size(newModel.S,2))]; +0335 +0336 if copyInfo==true +0337 [I, J]=ismember(metsToAdd.metNames,model.metNames); +0338 J=J(I); +0339 %I is the indexes of the new metabolites for which a metabolite with +0340 %the same name existed +0341 I=find(I)+nOldMets; +0342 %Go through each of the added mets and copy annotation if it doesn't +0343 %exist +0344 for i=1:numel(I) +0345 if isfield(newModel,'inchis') +0346 if isempty(newModel.inchis{I(i)}) +0347 newModel.inchis(I(i))=newModel.inchis(J(i)); +0348 end +0349 end +0350 if isfield(newModel,'metSmiles') +0351 if isempty(newModel.metSmiles{I(i)}) +0352 newModel.metSmiles(I(i))=newModel.metSmiles(J(i)); +0353 end +0354 end +0355 if isfield(newModel,'metFormulas') +0356 if isempty(newModel.metFormulas{I(i)}) +0357 newModel.metFormulas(I(i))=newModel.metFormulas(J(i)); +0358 end +0359 end +0360 if isfield(newModel,'metMiriams') +0361 if isempty(newModel.metMiriams{I(i)}) +0362 newModel.metMiriams(I(i))=newModel.metMiriams(J(i)); +0363 end +0364 end +0365 if isfield(newModel,'metCharges') +0366 newModel.metCharges(I(i))=newModel.metCharges(J(i)); +0367 end +0368 if isfield(newModel,'metDeltaG') +0369 newModel.metDeltaG(I(i))=newModel.metDeltaG(J(i)); +0370 end +0371 end +0372 end +0373 end
Generated by m2html © 2005
\ No newline at end of file From a999800ad5ebcca10a609a5054e6054c115dcea0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 28 Oct 2024 21:27:58 +0100 Subject: [PATCH 2/9] feat: changeGrRules multiple rxns same grRule --- core/changeGrRules.m | 3 ++ doc/core/changeGrRules.html | 59 +++++++++++++++++++------------------ 2 files changed, 34 insertions(+), 28 deletions(-) diff --git a/core/changeGrRules.m b/core/changeGrRules.m index 9933cb81..2ab838d3 100755 --- a/core/changeGrRules.m +++ b/core/changeGrRules.m @@ -23,6 +23,9 @@ rxns=convertCharArray(rxns); grRules=convertCharArray(grRules); +if isscalar(grRules) && ~isscalar(rxns) + grRules = repmat(grRules,1,numel(rxns)); +end if ~(numel(grRules)==numel(rxns)) error('Number of rxns and grRules should be identical') end diff --git a/doc/core/changeGrRules.html b/doc/core/changeGrRules.html index 2152eca1..95f7d385 100644 --- a/doc/core/changeGrRules.html +++ b/doc/core/changeGrRules.html @@ -82,35 +82,38 @@

SOURCE CODE ^convertCharArray(rxns); 0024 grRules=convertCharArray(grRules); 0025 -0026 if ~(numel(grRules)==numel(rxns)) -0027 error('Number of rxns and grRules should be identical') +0026 if isscalar(grRules) && ~isscalar(rxns) +0027 grRules = repmat(grRules,1,numel(rxns)); 0028 end -0029 -0030 for i=1:length(rxns) -0031 % Add genes to model -0032 geneList=transpose(cell(unique(regexp(grRules{i},'[)(]*|( and )*|( or )*','split')))); % Extract individual, unique genes from the geneAssoc provided -0033 geneList=geneList(~cellfun(@isempty, geneList)); -0034 genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes. -0035 if ~isempty(genesToAdd.genes) -0036 model=addGenesRaven(model,genesToAdd); % Add genes -0037 end -0038 -0039 % Find reaction and gene indices -0040 idx=getIndexes(model,rxns,'rxns'); -0041 end -0042 -0043 % Change gene associations -0044 if replace==true % Replace old gene associations -0045 model.grRules(idx)=grRules; -0046 else % Add gene associations, add new gene rules after 'OR'. -0047 model.grRules(idx)=strcat('(',model.grRules(idx),') or (',grRules,')'); -0048 end -0049 -0050 %Fix grRules and reconstruct rxnGeneMat -0051 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0052 model.grRules = grRules; -0053 model.rxnGeneMat = rxnGeneMat; -0054 end +0029 if ~(numel(grRules)==numel(rxns)) +0030 error('Number of rxns and grRules should be identical') +0031 end +0032 +0033 for i=1:length(rxns) +0034 % Add genes to model +0035 geneList=transpose(cell(unique(regexp(grRules{i},'[)(]*|( and )*|( or )*','split')))); % Extract individual, unique genes from the geneAssoc provided +0036 geneList=geneList(~cellfun(@isempty, geneList)); +0037 genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes. +0038 if ~isempty(genesToAdd.genes) +0039 model=addGenesRaven(model,genesToAdd); % Add genes +0040 end +0041 +0042 % Find reaction and gene indices +0043 idx=getIndexes(model,rxns,'rxns'); +0044 end +0045 +0046 % Change gene associations +0047 if replace==true % Replace old gene associations +0048 model.grRules(idx)=grRules; +0049 else % Add gene associations, add new gene rules after 'OR'. +0050 model.grRules(idx)=strcat('(',model.grRules(idx),') or (',grRules,')'); +0051 end +0052 +0053 %Fix grRules and reconstruct rxnGeneMat +0054 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0055 model.grRules = grRules; +0056 model.rxnGeneMat = rxnGeneMat; +0057 end
Generated by m2html © 2005
\ No newline at end of file From 20641923b770b91496a2a330a84a97f7e97898de Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 28 Oct 2024 21:28:39 +0100 Subject: [PATCH 3/9] fix: addRxns subSystems correct format --- core/addRxns.m | 2 +- doc/core/addRxns.html | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/core/addRxns.m b/core/addRxns.m index 1f7cf911..f23222e1 100755 --- a/core/addRxns.m +++ b/core/addRxns.m @@ -363,7 +363,7 @@ if ~isfield(newModel,'subSystems') newModel.subSystems=celllargefiller; end - newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems]; + newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)]; else %Fill with standard if it doesn't exist if isfield(newModel,'subSystems') diff --git a/doc/core/addRxns.html b/doc/core/addRxns.html index 11a2ec9b..3d0fc7a0 100644 --- a/doc/core/addRxns.html +++ b/doc/core/addRxns.html @@ -507,7 +507,7 @@

SOURCE CODE ^if ~isfield(newModel,'subSystems') 0364 newModel.subSystems=celllargefiller; 0365 end -0366 newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems]; +0366 newModel.subSystems=[newModel.subSystems;rxnsToAdd.subSystems(:)]; 0367 else 0368 %Fill with standard if it doesn't exist 0369 if isfield(newModel,'subSystems') From c90390e95cbb73ad75996f2e8b23f97ec84cea86 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 28 Oct 2024 21:32:23 +0100 Subject: [PATCH 4/9] fix: subSystems if cell array of cell arrays --- core/addExchangeRxns.m | 6 +++- doc/core/addExchangeRxns.html | 64 +++++++++++++++++++---------------- 2 files changed, 39 insertions(+), 31 deletions(-) diff --git a/core/addExchangeRxns.m b/core/addExchangeRxns.m index b4b6cf6f..31df9b80 100755 --- a/core/addExchangeRxns.m +++ b/core/addExchangeRxns.m @@ -62,7 +62,11 @@ model.eccodes=[model.eccodes;filler]; end if isfield(model,'subSystems') - model.subSystems=[model.subSystems;filler]; + fillerSub = filler; + if iscell(model.subSystems(1,1)) + fillerSub = repmat({fillerSub},numel(J),1); + end + model.subSystems=[model.subSystems;fillerSub]; end if isfield(model,'grRules') model.grRules=[model.grRules;filler]; diff --git a/doc/core/addExchangeRxns.html b/doc/core/addExchangeRxns.html index bdc91c13..77045cdd 100644 --- a/doc/core/addExchangeRxns.html +++ b/doc/core/addExchangeRxns.html @@ -126,37 +126,41 @@

SOURCE CODE ^end 0064 if isfield(model,'subSystems') -0065 model.subSystems=[model.subSystems;filler]; -0066 end -0067 if isfield(model,'grRules') -0068 model.grRules=[model.grRules;filler]; -0069 end -0070 if isfield(model,'rxnFrom') -0071 model.rxnFrom=[model.rxnFrom;filler]; -0072 end -0073 if isfield(model,'rxnMiriams') -0074 model.rxnMiriams=[model.rxnMiriams;filler]; -0075 end -0076 if isfield(model,'rxnGeneMat') -0077 model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(J),numel(model.genes))]; -0078 end -0079 if isfield(model,'rxnComps') -0080 model.rxnComps=[model.rxnComps;ones(numel(J),1)]; -0081 fprintf('NOTE: The exchange reactions are assigned to the first compartment\n'); +0065 fillerSub = filler; +0066 if iscell(model.subSystems(1,1)) +0067 fillerSub = repmat({fillerSub},numel(J),1); +0068 end +0069 model.subSystems=[model.subSystems;fillerSub]; +0070 end +0071 if isfield(model,'grRules') +0072 model.grRules=[model.grRules;filler]; +0073 end +0074 if isfield(model,'rxnFrom') +0075 model.rxnFrom=[model.rxnFrom;filler]; +0076 end +0077 if isfield(model,'rxnMiriams') +0078 model.rxnMiriams=[model.rxnMiriams;filler]; +0079 end +0080 if isfield(model,'rxnGeneMat') +0081 model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(J),numel(model.genes))]; 0082 end -0083 if isfield(model,'rxnNotes') -0084 model.rxnNotes=[model.rxnNotes;filler]; -0085 end -0086 if isfield(model,'rxnReferences') -0087 model.rxnReferences=[model.rxnReferences;filler]; -0088 end -0089 if isfield(model,'rxnConfidenceScores') -0090 model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(J),1)]; -0091 end -0092 if isfield(model,'rxnDeltaG') -0093 model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(J),1)]; -0094 end -0095 end +0083 if isfield(model,'rxnComps') +0084 model.rxnComps=[model.rxnComps;ones(numel(J),1)]; +0085 fprintf('NOTE: The exchange reactions are assigned to the first compartment\n'); +0086 end +0087 if isfield(model,'rxnNotes') +0088 model.rxnNotes=[model.rxnNotes;filler]; +0089 end +0090 if isfield(model,'rxnReferences') +0091 model.rxnReferences=[model.rxnReferences;filler]; +0092 end +0093 if isfield(model,'rxnConfidenceScores') +0094 model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(J),1)]; +0095 end +0096 if isfield(model,'rxnDeltaG') +0097 model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(J),1)]; +0098 end +0099 end
Generated by m2html © 2005
\ No newline at end of file From 3fbd2429e5dc3a29c9cbf81530bb83a0421c9814 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 29 Oct 2024 14:11:37 +0100 Subject: [PATCH 5/9] fix: importModel if identifiers.org/name:value --- doc/io/importModel.html | 91 +++++++++++++++++++++-------------------- io/importModel.m | 13 +++--- 2 files changed, 53 insertions(+), 51 deletions(-) diff --git a/doc/io/importModel.html b/doc/io/importModel.html index d0d442fc..ed93b7fc 100644 --- a/doc/io/importModel.html +++ b/doc/io/importModel.html @@ -1238,52 +1238,53 @@

SOURCE CODE ^%Finding whether miriams are written in the old or the new way 1129 if strfind(searchString,'urn:miriam:') 1130 startString='urn:miriam:'; -1131 midString=':'; -1132 elseif strfind(searchString,'http://identifiers.org/') -1133 startString='http://identifiers.org/'; -1134 midString='/'; -1135 elseif strfind(searchString,'https://identifiers.org/') -1136 startString='https://identifiers.org/'; -1137 midString='/'; -1138 else -1139 miriamStruct=[]; -1140 return; -1141 end -1142 -1143 miriamStruct=[]; -1144 -1145 searchString=regexprep(searchString,'" />','"/>'); -1146 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); -1147 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1148 targetString=regexprep(targetString,startString,''); -1149 targetString=regexprep(targetString,midString,'/','once'); -1150 -1151 counter=0; -1152 for i=1:numel(targetString) -1153 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) -1154 counter=counter+1; -1155 miriamStruct.name{counter,1} = regexprep(targetString{1,i},'/.+','','once'); -1156 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} '/'],'','once'); -1157 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); -1158 end -1159 end +1131 elseif strfind(searchString,'http://identifiers.org/') +1132 startString='http://identifiers.org/'; +1133 elseif strfind(searchString,'https://identifiers.org/') +1134 startString='https://identifiers.org/'; +1135 else +1136 miriamStruct=[]; +1137 return; +1138 end +1139 +1140 miriamStruct=[]; +1141 +1142 searchString=regexprep(searchString,'" />','"/>'); +1143 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); +1144 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1145 targetString=regexprep(targetString,startString,''); +1146 +1147 fwdslash = contains(targetString,'/'); +1148 midString = cell(numel(targetString),1); +1149 midString(fwdslash) = {'/'}; +1150 midString(~fwdslash) = {':'}; +1151 +1152 counter=0; +1153 for i=1:numel(targetString) +1154 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) +1155 counter=counter+1; +1156 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); +1157 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); +1158 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); +1159 end 1160 end -1161 -1162 function miriam = addSBOtoMiriam(miriam,sboTerm) -1163 %Appends SBO term to miriam structure -1164 -1165 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format -1166 if isempty(miriam) -1167 miriam.name = {'sbo'}; -1168 miriam.value = sboTerm; -1169 elseif any(strcmp('sbo',miriam.name)) -1170 currSbo = strcmp('sbo',miriam.name); -1171 miriam.value(currSbo) = sboTerm; -1172 else -1173 miriam.name(end+1) = {'sbo'}; -1174 miriam.value(end+1) = sboTerm; -1175 end -1176 end +1161 end +1162 +1163 function miriam = addSBOtoMiriam(miriam,sboTerm) +1164 %Appends SBO term to miriam structure +1165 +1166 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format +1167 if isempty(miriam) +1168 miriam.name = {'sbo'}; +1169 miriam.value = sboTerm; +1170 elseif any(strcmp('sbo',miriam.name)) +1171 currSbo = strcmp('sbo',miriam.name); +1172 miriam.value(currSbo) = sboTerm; +1173 else +1174 miriam.name(end+1) = {'sbo'}; +1175 miriam.value(end+1) = sboTerm; +1176 end +1177 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/io/importModel.m b/io/importModel.m index b36a4ad9..3032ba53 100755 --- a/io/importModel.m +++ b/io/importModel.m @@ -1128,13 +1128,10 @@ %Finding whether miriams are written in the old or the new way if strfind(searchString,'urn:miriam:') startString='urn:miriam:'; - midString=':'; elseif strfind(searchString,'http://identifiers.org/') startString='http://identifiers.org/'; - midString='/'; elseif strfind(searchString,'https://identifiers.org/') startString='https://identifiers.org/'; - midString='/'; else miriamStruct=[]; return; @@ -1146,14 +1143,18 @@ [~,targetString] = regexp(searchString,'','tokens','match'); targetString=regexprep(targetString,'',''); targetString=regexprep(targetString,startString,''); -targetString=regexprep(targetString,midString,'/','once'); + +fwdslash = contains(targetString,'/'); +midString = cell(numel(targetString),1); +midString(fwdslash) = {'/'}; +midString(~fwdslash) = {':'}; counter=0; for i=1:numel(targetString) if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) counter=counter+1; - miriamStruct.name{counter,1} = regexprep(targetString{1,i},'/.+','','once'); - miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} '/'],'','once'); + miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); + miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); end end From a9b1817de1b749cb6367ecdd362365814eb98dca Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 29 Oct 2024 14:12:28 +0100 Subject: [PATCH 6/9] refactor: checkModelStruct use getGenesFromGrRules --- core/checkModelStruct.m | 4 +- doc/core/checkModelStruct.html | 654 +++++++++++++++--------------- doc/core/getGenesFromGrRules.html | 2 +- 3 files changed, 328 insertions(+), 332 deletions(-) diff --git a/core/checkModelStruct.m b/core/checkModelStruct.m index 139c43d8..6cadafcf 100755 --- a/core/checkModelStruct.m +++ b/core/checkModelStruct.m @@ -124,9 +124,7 @@ function checkModelStruct(model,throwErrors,trimWarnings) EM='If "grRules" field exists, the model should also contain a "genes" field'; dispEM(EM,throwErrors); else - geneList = strjoin(model.grRules); - geneList = regexp(geneList,' |)|(|and|or','split'); % Remove all grRule punctuation - geneList = geneList(~cellfun(@isempty,geneList)); % Remove spaces and empty genes + geneList = getGenesFromGrRules(model.grRules); geneList = setdiff(unique(geneList),model.genes); if ~isempty(geneList) problemGrRules = model.rxns(contains(model.grRules,geneList)); diff --git a/doc/core/checkModelStruct.html b/doc/core/checkModelStruct.html index 288a1ff0..2f029e1f 100644 --- a/doc/core/checkModelStruct.html +++ b/doc/core/checkModelStruct.html @@ -46,7 +46,7 @@

DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

This function calls: +
  • dispEM dispEM
  • getGenesFromGrRules getGenesFromGrRules Extract gene list and rxnGeneMat from grRules array.
  • parseFormulas parseFormulas
    • This function is called by:
    @@ -183,333 +183,331 @@

    SOURCE CODE ^'If "grRules" field exists, the model should also contain a "genes" field'; 0125 dispEM(EM,throwErrors); 0126 else -0127 geneList = strjoin(model.grRules); -0128 geneList = regexp(geneList,' |)|(|and|or','split'); % Remove all grRule punctuation -0129 geneList = geneList(~cellfun(@isempty,geneList)); % Remove spaces and empty genes -0130 geneList = setdiff(unique(geneList),model.genes); -0131 if ~isempty(geneList) -0132 problemGrRules = model.rxns(contains(model.grRules,geneList)); -0133 problemGrRules = strjoin(problemGrRules(:),'; '); -0134 EM=['The reaction(s) "' problemGrRules '" contain the following genes in its "grRules" field, but these are not in the "genes" field:']; -0135 dispEM(EM,throwErrors,geneList); -0136 end -0137 end -0138 end -0139 if isfield(model,'rxnComps') -0140 if ~isnumeric(model.rxnComps) -0141 EM='The "rxnComps" field must be of type "double"'; -0142 dispEM(EM,throwErrors); -0143 end -0144 end -0145 if isfield(model,'inchis') -0146 if ~iscellstr(model.inchis) -0147 EM='The "inchis" field must be a cell array of strings'; -0148 dispEM(EM,throwErrors); -0149 end -0150 end -0151 if isfield(model,'metSmiles') -0152 if ~iscellstr(model.metSmiles) -0153 EM='The "metSmiles" field must be a cell array of strings'; -0154 dispEM(EM,throwErrors); -0155 end -0156 end -0157 if isfield(model,'metFormulas') -0158 if ~iscellstr(model.metFormulas) -0159 EM='The "metFormulas" field must be a cell array of strings'; -0160 dispEM(EM,throwErrors); -0161 end -0162 end -0163 if isfield(model,'metCharges') -0164 if ~isnumeric(model.metCharges) -0165 EM='The "metCharges" field must be a double'; -0166 dispEM(EM,throwErrors); -0167 end -0168 end -0169 if isfield(model,'metDeltaG') -0170 if ~isnumeric(model.metDeltaG) -0171 EM='The "metDeltaG" field must be a double'; -0172 dispEM(EM,throwErrors); -0173 end -0174 end -0175 if isfield(model,'subSystems') -0176 for i=1:numel(model.subSystems) -0177 if ~iscell(model.subSystems{i,1}) -0178 EM='The "subSystems" field must be a cell array'; -0179 dispEM(EM,throwErrors); -0180 end -0181 end -0182 end -0183 if isfield(model,'eccodes') -0184 if ~iscellstr(model.eccodes) -0185 EM='The "eccodes" field must be a cell array of strings'; -0186 dispEM(EM,throwErrors); -0187 end -0188 end -0189 if isfield(model,'unconstrained') -0190 if ~isnumeric(model.unconstrained) -0191 EM='The "unconstrained" field must be of type "double"'; -0192 dispEM(EM,throwErrors); -0193 end -0194 end -0195 if isfield(model,'rxnNotes') -0196 if ~iscellstr(model.rxnNotes) -0197 EM='The "rxnNotes" field must be a cell array of strings'; -0198 dispEM(EM,throwErrors); -0199 end -0200 end -0201 if isfield(model,'rxnReferences') -0202 if ~iscellstr(model.rxnReferences) -0203 EM='The "rxnReferences" field must be a cell array of strings'; -0204 dispEM(EM,throwErrors); -0205 end -0206 end -0207 if isfield(model,'rxnConfidenceScores') -0208 if ~isnumeric(model.rxnConfidenceScores) -0209 EM='The "rxnConfidenceScores" field must be a double'; -0210 dispEM(EM,throwErrors); -0211 end -0212 end -0213 if isfield(model,'rxnDeltaG') -0214 if ~isnumeric(model.rxnDeltaG) -0215 EM='The "rxnDeltaG" field must be a double'; -0216 dispEM(EM,throwErrors); -0217 end -0218 end -0219 -0220 %Empty strings -0221 if isempty(model.id) -0222 EM='The "id" field cannot be empty'; -0223 dispEM(EM,throwErrors); -0224 end -0225 if any(cellfun(@isempty,model.rxns)) -0226 EM='The model contains empty reaction IDs'; -0227 dispEM(EM,throwErrors); -0228 end -0229 if any(cellfun(@isempty,model.mets)) -0230 EM='The model contains empty metabolite IDs'; -0231 dispEM(EM,throwErrors); -0232 end -0233 if any(cellfun(@isempty,model.comps)) -0234 EM='The model contains empty compartment IDs'; -0235 dispEM(EM,throwErrors); -0236 end -0237 EM='The following metabolites have empty names:'; -0238 dispEM(EM,throwErrors,model.mets(cellfun(@isempty,model.metNames)),trimWarnings); -0239 -0240 if isfield(model,'genes') -0241 if any(cellfun(@isempty,model.genes)) -0242 EM='The model contains empty gene IDs'; -0243 dispEM(EM,throwErrors); -0244 end -0245 end -0246 -0247 %Validate format of ids -0248 fields = {'rxns';'mets';'comps';'genes'}; -0249 fieldNames = {'reaction';'metabolite';'compartment';'gene'}; -0250 fieldPrefix = {'R_';'M_';'C_';'G_'}; -0251 for i=1:numel(fields) -0252 try -0253 numIDs = ~startsWith(model.(fields{i}),regexpPattern('^[a-zA-Z_]')); -0254 catch -0255 numIDs = []; -0256 end -0257 if any(numIDs) -0258 EM = ['The following ' fieldNames{i} ' identifiers do not start '... -0259 'with a letter or _ (conflicting with SBML specifications). '... -0260 'This does not impact RAVEN functionality, but be aware that '... -0261 'exportModel will automatically add ' fieldPrefix{i} ... -0262 ' prefixes to all ' fieldNames{i} ' identifiers:']; -0263 dispEM(EM,false,{model.(fields{i}){numIDs}},trimWarnings); -0264 end -0265 end -0266 -0267 %Duplicates -0268 EM='The following reaction IDs are duplicates:'; -0269 dispEM(EM,throwErrors,model.rxns(duplicates(model.rxns)),trimWarnings); -0270 EM='The following metabolite IDs are duplicates:'; -0271 dispEM(EM,throwErrors,model.mets(duplicates(model.mets)),trimWarnings); -0272 EM='The following compartment IDs are duplicates:'; -0273 dispEM(EM,throwErrors,model.comps(duplicates(model.comps)),trimWarnings); -0274 if isfield(model,'genes') -0275 EM='The following genes are duplicates:'; -0276 dispEM(EM,throwErrors,model.genes(duplicates(model.genes)),trimWarnings); -0277 end -0278 metInComp=strcat(model.metNames,'[',model.comps(model.metComps),']'); -0279 EM='The following metabolites already exist in the same compartment:'; -0280 dispEM(EM,throwErrors,metInComp(duplicates(metInComp)),trimWarnings); -0281 -0282 %Elements never used (print only as warnings -0283 EM='The following reactions are empty (no involved metabolites):'; -0284 dispEM(EM,false,model.rxns(~any(model.S,1)),trimWarnings); -0285 EM='The following metabolites are never used in a reaction:'; -0286 dispEM(EM,false,model.mets(~any(model.S,2)),trimWarnings); -0287 if isfield(model,'genes') -0288 EM='The following genes are not associated to a reaction:'; -0289 dispEM(EM,false,model.genes(~any(model.rxnGeneMat,1)),trimWarnings); -0290 end -0291 I=true(numel(model.comps),1); -0292 I(model.metComps)=false; -0293 EM='The following compartments contain no metabolites:'; -0294 dispEM(EM,false,model.comps(I),trimWarnings); -0295 -0296 %Contradicting bounds -0297 EM='The following reactions have contradicting bounds (lower bound is higher than upper bound):'; -0298 dispEM(EM,throwErrors,model.rxns(model.lb>model.ub),trimWarnings); -0299 EM='The following reactions have lower and upper bounds that indicate reversibility, but are indicated as irreversible in model.rev:'; -0300 dispEM(EM,false,model.rxns(model.lb < 0 & model.ub > 0 & model.rev==0),trimWarnings); -0301 -0302 %Multiple or no objective functions not allowed in SBML L3V1 FBCv2 -0303 if numel(find(model.c))>1 -0304 EM='Multiple objective functions found. This might be intended, but results in FBCv2 non-compliant SBML file when exported'; -0305 dispEM(EM,false,model.rxns(find(model.c)),trimWarnings); -0306 elseif ~any(model.c) -0307 EM='No objective function found. This might be intended, but results in FBCv2 non-compliant SBML file when exported'; -0308 dispEM(EM,false); -0309 end -0310 -0311 %Mapping of compartments -0312 if isfield(model,'compOutside') -0313 EM='The following compartments are in "compOutside" but not in "comps":'; -0314 dispEM(EM,throwErrors,setdiff(model.compOutside,[{''};model.comps]),trimWarnings); -0315 end -0316 -0317 %Met names which start with number -0318 I=false(numel(model.metNames),1); -0319 for i=1:numel(model.metNames) -0320 index=strfind(model.metNames{i},' '); -0321 if any(index) -0322 if any(str2double(model.metNames{i}(1:index(1)-1))) -0323 I(i)=true; -0324 end -0325 end -0326 end -0327 EM='The following metabolite names begin with a number directly followed by space, which could potentially cause problems:'; -0328 dispEM(EM,false,model.metNames(I),trimWarnings); -0329 -0330 %Non-parseable composition -0331 if isfield(model,'metFormulas') -0332 [~, ~, exitFlag]=parseFormulas(model.metFormulas,true,false); -0333 EM='The composition for the following metabolites could not be parsed:'; -0334 dispEM(EM,false,model.mets(exitFlag==-1),trimWarnings); -0335 end -0336 -0337 %Check if there are metabolites with different names but the same MIRIAM -0338 %codes -0339 if isfield(model,'metMiriams') -0340 miriams=containers.Map(); -0341 for i=1:numel(model.mets) -0342 if ~isempty(model.metMiriams{i}) -0343 %Loop through and add for each miriam -0344 for j=1:numel(model.metMiriams{i}.name) -0345 %Get existing metabolite indexes -0346 current=strcat(model.metMiriams{i}.name{j},'/',model.metMiriams{i}.value{j}); -0347 if isKey(miriams,current) -0348 existing=miriams(current); -0349 else -0350 existing=[]; -0351 end -0352 miriams(current)=[existing;i]; -0353 end -0354 end -0355 end -0356 -0357 %Get all keys -0358 allMiriams=keys(miriams); -0359 -0360 hasMultiple=false(numel(allMiriams),1); -0361 for i=1:numel(allMiriams) -0362 if numel(miriams(allMiriams{i}))>1 -0363 %Check if they all have the same name -0364 if numel(unique(model.metNames(miriams(allMiriams{i}))))>1 -0365 if ~regexp(allMiriams{i},'^sbo\/SBO:') % SBO terms are expected to be multiple -0366 hasMultiple(i)=true; -0367 end -0368 end -0369 end -0370 end -0371 -0372 %Print output -0373 EM='The following MIRIAM strings are associated to more than one unique metabolite name:'; -0374 dispEM(EM,false,allMiriams(hasMultiple),trimWarnings); -0375 end -0376 -0377 %Check if there are metabolites with different names but the same InChI -0378 %codes -0379 if isfield(model,'inchis') -0380 inchis=containers.Map(); -0381 for i=1:numel(model.mets) -0382 if ~isempty(model.inchis{i}) -0383 %Get existing metabolite indexes -0384 if isKey(inchis,model.inchis{i}) -0385 existing=inchis(model.inchis{i}); -0386 else -0387 existing=[]; -0388 end -0389 inchis(model.inchis{i})=[existing;i]; -0390 end -0391 end -0392 -0393 %Get all keys -0394 allInchis=keys(inchis); -0395 -0396 hasMultiple=false(numel(allInchis),1); -0397 for i=1:numel(allInchis) -0398 if numel(inchis(allInchis{i}))>1 -0399 %Check if they all have the same name -0400 if numel(unique(model.metNames(inchis(allInchis{i}))))>1 -0401 hasMultiple(i)=true; -0402 end -0403 end -0404 end -0405 -0406 %Print output -0407 EM='The following InChI strings are associated to more than one unique metabolite name:'; -0408 dispEM(EM,false,allInchis(hasMultiple),trimWarnings); -0409 end -0410 -0411 % %Check if there are metabolites with different names but the same SMILES -0412 % if isfield(model,'metSmiles') -0413 % metSmiles=containers.Map(); -0414 % for i=1:numel(model.mets) -0415 % if ~isempty(model.metSmiles{i}) -0416 % %Get existing metabolite indexes -0417 % if isKey(metSmiles,model.metSmiles{i}) -0418 % existing=metSmiles(model.metSmiles{i}); -0419 % else -0420 % existing=[]; -0421 % end -0422 % metSmiles(model.metSmiles{i})=[existing;i]; -0423 % end -0424 % end -0425 % -0426 % %Get all keys -0427 % allmetSmiles=keys(metSmiles); -0428 % -0429 % hasMultiple=false(numel(metSmiles),1); -0430 % for i=1:numel(metSmiles) -0431 % if numel(metSmiles(metSmiles{i}))>1 -0432 % %Check if they all have the same name -0433 % if numel(unique(model.metNames(metSmiles(allmetSmiles{i}))))>1 -0434 % hasMultiple(i)=true; -0435 % end -0436 % end -0437 % end -0438 % -0439 % %Print output -0440 % EM='The following metSmiles strings are associated to more than one unique metabolite name:'; -0441 % dispEM(EM,false,allmetSmiles(hasMultiple),trimWarnings); -0442 % end -0443 end -0444 -0445 function I=duplicates(strings) -0446 I=false(numel(strings),1); -0447 [J, K]=unique(strings); -0448 if numel(J)~=numel(strings) -0449 L=1:numel(strings); -0450 L(K)=[]; -0451 I(L)=true; -0452 end -0453 end +0127 geneList = getGenesFromGrRules(model.grRules); +0128 geneList = setdiff(unique(geneList),model.genes); +0129 if ~isempty(geneList) +0130 problemGrRules = model.rxns(contains(model.grRules,geneList)); +0131 problemGrRules = strjoin(problemGrRules(:),'; '); +0132 EM=['The reaction(s) "' problemGrRules '" contain the following genes in its "grRules" field, but these are not in the "genes" field:']; +0133 dispEM(EM,throwErrors,geneList); +0134 end +0135 end +0136 end +0137 if isfield(model,'rxnComps') +0138 if ~isnumeric(model.rxnComps) +0139 EM='The "rxnComps" field must be of type "double"'; +0140 dispEM(EM,throwErrors); +0141 end +0142 end +0143 if isfield(model,'inchis') +0144 if ~iscellstr(model.inchis) +0145 EM='The "inchis" field must be a cell array of strings'; +0146 dispEM(EM,throwErrors); +0147 end +0148 end +0149 if isfield(model,'metSmiles') +0150 if ~iscellstr(model.metSmiles) +0151 EM='The "metSmiles" field must be a cell array of strings'; +0152 dispEM(EM,throwErrors); +0153 end +0154 end +0155 if isfield(model,'metFormulas') +0156 if ~iscellstr(model.metFormulas) +0157 EM='The "metFormulas" field must be a cell array of strings'; +0158 dispEM(EM,throwErrors); +0159 end +0160 end +0161 if isfield(model,'metCharges') +0162 if ~isnumeric(model.metCharges) +0163 EM='The "metCharges" field must be a double'; +0164 dispEM(EM,throwErrors); +0165 end +0166 end +0167 if isfield(model,'metDeltaG') +0168 if ~isnumeric(model.metDeltaG) +0169 EM='The "metDeltaG" field must be a double'; +0170 dispEM(EM,throwErrors); +0171 end +0172 end +0173 if isfield(model,'subSystems') +0174 for i=1:numel(model.subSystems) +0175 if ~iscell(model.subSystems{i,1}) +0176 EM='The "subSystems" field must be a cell array'; +0177 dispEM(EM,throwErrors); +0178 end +0179 end +0180 end +0181 if isfield(model,'eccodes') +0182 if ~iscellstr(model.eccodes) +0183 EM='The "eccodes" field must be a cell array of strings'; +0184 dispEM(EM,throwErrors); +0185 end +0186 end +0187 if isfield(model,'unconstrained') +0188 if ~isnumeric(model.unconstrained) +0189 EM='The "unconstrained" field must be of type "double"'; +0190 dispEM(EM,throwErrors); +0191 end +0192 end +0193 if isfield(model,'rxnNotes') +0194 if ~iscellstr(model.rxnNotes) +0195 EM='The "rxnNotes" field must be a cell array of strings'; +0196 dispEM(EM,throwErrors); +0197 end +0198 end +0199 if isfield(model,'rxnReferences') +0200 if ~iscellstr(model.rxnReferences) +0201 EM='The "rxnReferences" field must be a cell array of strings'; +0202 dispEM(EM,throwErrors); +0203 end +0204 end +0205 if isfield(model,'rxnConfidenceScores') +0206 if ~isnumeric(model.rxnConfidenceScores) +0207 EM='The "rxnConfidenceScores" field must be a double'; +0208 dispEM(EM,throwErrors); +0209 end +0210 end +0211 if isfield(model,'rxnDeltaG') +0212 if ~isnumeric(model.rxnDeltaG) +0213 EM='The "rxnDeltaG" field must be a double'; +0214 dispEM(EM,throwErrors); +0215 end +0216 end +0217 +0218 %Empty strings +0219 if isempty(model.id) +0220 EM='The "id" field cannot be empty'; +0221 dispEM(EM,throwErrors); +0222 end +0223 if any(cellfun(@isempty,model.rxns)) +0224 EM='The model contains empty reaction IDs'; +0225 dispEM(EM,throwErrors); +0226 end +0227 if any(cellfun(@isempty,model.mets)) +0228 EM='The model contains empty metabolite IDs'; +0229 dispEM(EM,throwErrors); +0230 end +0231 if any(cellfun(@isempty,model.comps)) +0232 EM='The model contains empty compartment IDs'; +0233 dispEM(EM,throwErrors); +0234 end +0235 EM='The following metabolites have empty names:'; +0236 dispEM(EM,throwErrors,model.mets(cellfun(@isempty,model.metNames)),trimWarnings); +0237 +0238 if isfield(model,'genes') +0239 if any(cellfun(@isempty,model.genes)) +0240 EM='The model contains empty gene IDs'; +0241 dispEM(EM,throwErrors); +0242 end +0243 end +0244 +0245 %Validate format of ids +0246 fields = {'rxns';'mets';'comps';'genes'}; +0247 fieldNames = {'reaction';'metabolite';'compartment';'gene'}; +0248 fieldPrefix = {'R_';'M_';'C_';'G_'}; +0249 for i=1:numel(fields) +0250 try +0251 numIDs = ~startsWith(model.(fields{i}),regexpPattern('^[a-zA-Z_]')); +0252 catch +0253 numIDs = []; +0254 end +0255 if any(numIDs) +0256 EM = ['The following ' fieldNames{i} ' identifiers do not start '... +0257 'with a letter or _ (conflicting with SBML specifications). '... +0258 'This does not impact RAVEN functionality, but be aware that '... +0259 'exportModel will automatically add ' fieldPrefix{i} ... +0260 ' prefixes to all ' fieldNames{i} ' identifiers:']; +0261 dispEM(EM,false,{model.(fields{i}){numIDs}},trimWarnings); +0262 end +0263 end +0264 +0265 %Duplicates +0266 EM='The following reaction IDs are duplicates:'; +0267 dispEM(EM,throwErrors,model.rxns(duplicates(model.rxns)),trimWarnings); +0268 EM='The following metabolite IDs are duplicates:'; +0269 dispEM(EM,throwErrors,model.mets(duplicates(model.mets)),trimWarnings); +0270 EM='The following compartment IDs are duplicates:'; +0271 dispEM(EM,throwErrors,model.comps(duplicates(model.comps)),trimWarnings); +0272 if isfield(model,'genes') +0273 EM='The following genes are duplicates:'; +0274 dispEM(EM,throwErrors,model.genes(duplicates(model.genes)),trimWarnings); +0275 end +0276 metInComp=strcat(model.metNames,'[',model.comps(model.metComps),']'); +0277 EM='The following metabolites already exist in the same compartment:'; +0278 dispEM(EM,throwErrors,metInComp(duplicates(metInComp)),trimWarnings); +0279 +0280 %Elements never used (print only as warnings +0281 EM='The following reactions are empty (no involved metabolites):'; +0282 dispEM(EM,false,model.rxns(~any(model.S,1)),trimWarnings); +0283 EM='The following metabolites are never used in a reaction:'; +0284 dispEM(EM,false,model.mets(~any(model.S,2)),trimWarnings); +0285 if isfield(model,'genes') +0286 EM='The following genes are not associated to a reaction:'; +0287 dispEM(EM,false,model.genes(~any(model.rxnGeneMat,1)),trimWarnings); +0288 end +0289 I=true(numel(model.comps),1); +0290 I(model.metComps)=false; +0291 EM='The following compartments contain no metabolites:'; +0292 dispEM(EM,false,model.comps(I),trimWarnings); +0293 +0294 %Contradicting bounds +0295 EM='The following reactions have contradicting bounds (lower bound is higher than upper bound):'; +0296 dispEM(EM,throwErrors,model.rxns(model.lb>model.ub),trimWarnings); +0297 EM='The following reactions have lower and upper bounds that indicate reversibility, but are indicated as irreversible in model.rev:'; +0298 dispEM(EM,false,model.rxns(model.lb < 0 & model.ub > 0 & model.rev==0),trimWarnings); +0299 +0300 %Multiple or no objective functions not allowed in SBML L3V1 FBCv2 +0301 if numel(find(model.c))>1 +0302 EM='Multiple objective functions found. This might be intended, but results in FBCv2 non-compliant SBML file when exported'; +0303 dispEM(EM,false,model.rxns(find(model.c)),trimWarnings); +0304 elseif ~any(model.c) +0305 EM='No objective function found. This might be intended, but results in FBCv2 non-compliant SBML file when exported'; +0306 dispEM(EM,false); +0307 end +0308 +0309 %Mapping of compartments +0310 if isfield(model,'compOutside') +0311 EM='The following compartments are in "compOutside" but not in "comps":'; +0312 dispEM(EM,throwErrors,setdiff(model.compOutside,[{''};model.comps]),trimWarnings); +0313 end +0314 +0315 %Met names which start with number +0316 I=false(numel(model.metNames),1); +0317 for i=1:numel(model.metNames) +0318 index=strfind(model.metNames{i},' '); +0319 if any(index) +0320 if any(str2double(model.metNames{i}(1:index(1)-1))) +0321 I(i)=true; +0322 end +0323 end +0324 end +0325 EM='The following metabolite names begin with a number directly followed by space, which could potentially cause problems:'; +0326 dispEM(EM,false,model.metNames(I),trimWarnings); +0327 +0328 %Non-parseable composition +0329 if isfield(model,'metFormulas') +0330 [~, ~, exitFlag]=parseFormulas(model.metFormulas,true,false); +0331 EM='The composition for the following metabolites could not be parsed:'; +0332 dispEM(EM,false,model.mets(exitFlag==-1),trimWarnings); +0333 end +0334 +0335 %Check if there are metabolites with different names but the same MIRIAM +0336 %codes +0337 if isfield(model,'metMiriams') +0338 miriams=containers.Map(); +0339 for i=1:numel(model.mets) +0340 if ~isempty(model.metMiriams{i}) +0341 %Loop through and add for each miriam +0342 for j=1:numel(model.metMiriams{i}.name) +0343 %Get existing metabolite indexes +0344 current=strcat(model.metMiriams{i}.name{j},'/',model.metMiriams{i}.value{j}); +0345 if isKey(miriams,current) +0346 existing=miriams(current); +0347 else +0348 existing=[]; +0349 end +0350 miriams(current)=[existing;i]; +0351 end +0352 end +0353 end +0354 +0355 %Get all keys +0356 allMiriams=keys(miriams); +0357 +0358 hasMultiple=false(numel(allMiriams),1); +0359 for i=1:numel(allMiriams) +0360 if numel(miriams(allMiriams{i}))>1 +0361 %Check if they all have the same name +0362 if numel(unique(model.metNames(miriams(allMiriams{i}))))>1 +0363 if ~regexp(allMiriams{i},'^sbo\/SBO:') % SBO terms are expected to be multiple +0364 hasMultiple(i)=true; +0365 end +0366 end +0367 end +0368 end +0369 +0370 %Print output +0371 EM='The following MIRIAM strings are associated to more than one unique metabolite name:'; +0372 dispEM(EM,false,allMiriams(hasMultiple),trimWarnings); +0373 end +0374 +0375 %Check if there are metabolites with different names but the same InChI +0376 %codes +0377 if isfield(model,'inchis') +0378 inchis=containers.Map(); +0379 for i=1:numel(model.mets) +0380 if ~isempty(model.inchis{i}) +0381 %Get existing metabolite indexes +0382 if isKey(inchis,model.inchis{i}) +0383 existing=inchis(model.inchis{i}); +0384 else +0385 existing=[]; +0386 end +0387 inchis(model.inchis{i})=[existing;i]; +0388 end +0389 end +0390 +0391 %Get all keys +0392 allInchis=keys(inchis); +0393 +0394 hasMultiple=false(numel(allInchis),1); +0395 for i=1:numel(allInchis) +0396 if numel(inchis(allInchis{i}))>1 +0397 %Check if they all have the same name +0398 if numel(unique(model.metNames(inchis(allInchis{i}))))>1 +0399 hasMultiple(i)=true; +0400 end +0401 end +0402 end +0403 +0404 %Print output +0405 EM='The following InChI strings are associated to more than one unique metabolite name:'; +0406 dispEM(EM,false,allInchis(hasMultiple),trimWarnings); +0407 end +0408 +0409 % %Check if there are metabolites with different names but the same SMILES +0410 % if isfield(model,'metSmiles') +0411 % metSmiles=containers.Map(); +0412 % for i=1:numel(model.mets) +0413 % if ~isempty(model.metSmiles{i}) +0414 % %Get existing metabolite indexes +0415 % if isKey(metSmiles,model.metSmiles{i}) +0416 % existing=metSmiles(model.metSmiles{i}); +0417 % else +0418 % existing=[]; +0419 % end +0420 % metSmiles(model.metSmiles{i})=[existing;i]; +0421 % end +0422 % end +0423 % +0424 % %Get all keys +0425 % allmetSmiles=keys(metSmiles); +0426 % +0427 % hasMultiple=false(numel(metSmiles),1); +0428 % for i=1:numel(metSmiles) +0429 % if numel(metSmiles(metSmiles{i}))>1 +0430 % %Check if they all have the same name +0431 % if numel(unique(model.metNames(metSmiles(allmetSmiles{i}))))>1 +0432 % hasMultiple(i)=true; +0433 % end +0434 % end +0435 % end +0436 % +0437 % %Print output +0438 % EM='The following metSmiles strings are associated to more than one unique metabolite name:'; +0439 % dispEM(EM,false,allmetSmiles(hasMultiple),trimWarnings); +0440 % end +0441 end +0442 +0443 function I=duplicates(strings) +0444 I=false(numel(strings),1); +0445 [J, K]=unique(strings); +0446 if numel(J)~=numel(strings) +0447 L=1:numel(strings); +0448 L(K)=[]; +0449 I(L)=true; +0450 end +0451 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/core/getGenesFromGrRules.html b/doc/core/getGenesFromGrRules.html index 84f78e8c..15f2a098 100644 --- a/doc/core/getGenesFromGrRules.html +++ b/doc/core/getGenesFromGrRules.html @@ -58,7 +58,7 @@

    CROSS-REFERENCE INFORMATION ^ </ul> This function is called by: <ul style= - +
  • checkModelStruct checkModelStruct
    • From b548aa59d56caabf1c61a16e52a33184efb4551a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 29 Oct 2024 14:14:29 +0100 Subject: [PATCH 7/9] feat: replaceMets can use identifiers instead --- core/replaceMets.m | 79 +++++++---- doc/core/contractModel.html | 2 +- doc/core/replaceMets.html | 253 ++++++++++++++++++++---------------- 3 files changed, 194 insertions(+), 140 deletions(-) diff --git a/core/replaceMets.m b/core/replaceMets.m index 5d0850f3..9301d862 100755 --- a/core/replaceMets.m +++ b/core/replaceMets.m @@ -1,4 +1,4 @@ -function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose) +function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose,identifiers) % replaceMets % Replaces metabolite names and annotation with replacement metabolite % that is already in the model. If this results in duplicate metabolites, @@ -13,6 +13,9 @@ % verbose logical whether to print the ids of reactions that % involve the replaced metabolite (optional, default % false) +% identifiers true if 'metabolite' and 'replacement' refer to +% metabolite identifiers instead of metabolite names +% (optional, default false) % % Output: % model model structure with selected metabolites replaced @@ -20,35 +23,49 @@ % idxDuplRxns index of removedRxns in original model % % Note: This function is useful when the model contains both 'oxygen' and -% 'o2' as metabolites. +% 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, +% then the 'identifiers' flag should be set to true. % % Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose) metabolite=char(metabolite); replacement=char(replacement); -if nargin<4 +if nargin<4 || isempty(verbose) verbose=false; end +if nargin<5 + identifiers = false; +end % Find occurence of replacement metabolites. Annotation will be taken from -% first metabolite found. Metabolite ID from replacement will be used where -% possible. -repIdx = find(strcmp(replacement,model.metNames)); +% first metabolite found. +if identifiers + repIdx = find(strcmp(replacement,model.mets)); +else + repIdx = find(strcmp(replacement,model.metNames)); +end if isempty(repIdx) - error('The replacement metabolite name cannot be found in the model.'); + error('The replacement metabolite cannot be found in the model.'); end % Change name and information from metabolite to replacement metabolite -metIdx = find(strcmp(metabolite,model.metNames)); +if identifiers + metIdx = find(strcmp(metabolite,model.mets)); +else + metIdx = find(strcmp(metabolite,model.metNames)); +end if isempty(metIdx) - error('The to-be-replaced metabolite name cannot be found in the model.'); + error('The to-be-replaced metabolite cannot be found in the model.'); end + +rxnsWithMet = find(model.S(metIdx,:)); if verbose==true - fprintf('\n\nThe following reactions contain the replaced metabolite as reactant:\n') - fprintf(strjoin(model.rxns(find(model.S(metIdx,:))),'\n')) + fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') + fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) fprintf('\n') end + model.metNames(metIdx) = model.metNames(repIdx(1)); if isfield(model,'metFormulas') model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); @@ -68,24 +85,32 @@ if isfield(model,'metSmiles') model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); end -% Run through replacement metabolites and their compartments. If any of the -% to-be-replaced metabolites is already present (checked by -% metaboliteName[compartment], then the replacement metabolite is kept and -% the to-be-replace metabolite ID deleted. - -% Build list of metaboliteName[compartment] -metCompsN =cellstr(num2str(model.metComps)); -map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); -metCompsN = map.values(metCompsN); -metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); idxDelete=[]; -for i = 1:length(repIdx) - metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); - if length(metCompsNidx)>1 - for j = 2:length(metCompsNidx) - model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); - idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete +if identifiers + originalStoch = model.S(metIdx,rxnsWithMet); + model.S(repIdx,rxnsWithMet) = originalStoch; + model.S(metIdx,rxnsWithMet) = 0; + idxDelete = metIdx; +else + % Run through replacement metabolites and their compartments. If any of the + % to-be-replaced metabolites is already present (checked by + % metaboliteName[compartment], then the replacement metabolite is kept and + % the to-be-replace metabolite ID deleted. + + % Build list of metaboliteName[compartment] + metCompsN =cellstr(num2str(model.metComps)); + map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); + metCompsN = map.values(metCompsN); + metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); + + for i = 1:length(repIdx) + metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); + if length(metCompsNidx)>1 + for j = 2:length(metCompsNidx) + model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); + idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete + end end end end diff --git a/doc/core/contractModel.html b/doc/core/contractModel.html index 00b0e934..260bdc3e 100644 --- a/doc/core/contractModel.html +++ b/doc/core/contractModel.html @@ -62,7 +62,7 @@

    CROSS-REFERENCE INFORMATION ^ <li><a href=dispEM dispEM
  • getIndexes getIndexes
  • removeReactions removeReactions
  • sortModel sortModel
  • standardizeGrRules standardizeGrRules
    • This function is called by: +
  • mergeCompartments mergeCompartments
  • replaceMets replaceMets
  • simplifyModel simplifyModel

    • SUBFUNCTIONS ^

    diff --git a/doc/core/replaceMets.html b/doc/core/replaceMets.html index 937c095c..6f431951 100644 --- a/doc/core/replaceMets.html +++ b/doc/core/replaceMets.html @@ -24,7 +24,7 @@

    PURPOSE ^replaceMets

    SYNOPSIS ^

    -
    function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose)
    +
    function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose,identifiers)

    DESCRIPTION ^

     replaceMets
@@ -41,6 +41,9 @@ 
    DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
 
 
 <h2><a name=SOURCE CODE ^
    
-
    0001 function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose)
+
    0001 function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose,identifiers)
 0002 % replaceMets
 0003 %   Replaces metabolite names and annotation with replacement metabolite
 0004 %   that is already in the model. If this results in duplicate metabolites,
@@ -80,118 +84,143 @@ 
    SOURCE CODE ^%   verbose         logical whether to print the ids of reactions that
 0014 %                   involve the replaced metabolite (optional, default
 0015 %                   false)
-0016 %
-0017 % Output:
-0018 %   model           model structure with selected metabolites replaced
-0019 %   removedRxns     identifiers of duplicate reactions that were removed
-0020 %   idxDuplRxns     index of removedRxns in original model
-0021 %
-0022 % Note: This function is useful when the model contains both 'oxygen' and
-0023 % 'o2' as metabolites.
+0016 %   identifiers     true if 'metabolite' and 'replacement' refer to
+0017 %                   metabolite identifiers instead of metabolite names
+0018 %                   (optional, default false)
+0019 %
+0020 % Output:
+0021 %   model           model structure with selected metabolites replaced
+0022 %   removedRxns     identifiers of duplicate reactions that were removed
+0023 %   idxDuplRxns     index of removedRxns in original model
 0024 %
-0025 % Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose)
-0026 
-0027 metabolite=char(metabolite);
-0028 replacement=char(replacement);
-0029 
-0030 if nargin<4
-0031     verbose=false;
-0032 end
+0025 % Note: This function is useful when the model contains both 'oxygen' and
+0026 % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead,
+0027 % then the 'identifiers' flag should be set to true.
+0028 %
+0029 % Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose)
+0030 
+0031 metabolite=char(metabolite);
+0032 replacement=char(replacement);
 0033 
-0034 % Find occurence of replacement metabolites. Annotation will be taken from
-0035 % first metabolite found. Metabolite ID from replacement will be used where
-0036 % possible.
-0037 repIdx = find(strcmp(replacement,model.metNames));
-0038 if isempty(repIdx)
-0039     error('The replacement metabolite name cannot be found in the model.');
-0040 end
-0041 
-0042 % Change name and information from metabolite to replacement metabolite
-0043 metIdx = find(strcmp(metabolite,model.metNames));
-0044 if isempty(metIdx)
-0045     error('The to-be-replaced metabolite name cannot be found in the model.');
-0046 end
-0047 if verbose==true
-0048     fprintf('\n\nThe following reactions contain the replaced metabolite as reactant:\n')
-0049     fprintf(strjoin(model.rxns(find(model.S(metIdx,:))),'\n'))
-0050     fprintf('\n')
-0051 end
-0052 model.metNames(metIdx) = model.metNames(repIdx(1));
-0053 if isfield(model,'metFormulas')
-0054     model.metFormulas(metIdx) = model.metFormulas(repIdx(1));
-0055 end
-0056 if isfield(model,'metMiriams')
-0057     model.metMiriams(metIdx) = model.metMiriams(repIdx(1));
-0058 end
-0059 if isfield(model,'metCharges')
-0060     model.metCharges(metIdx) = model.metCharges(repIdx(1));
-0061 end
-0062 if isfield(model,'metDeltaG')
-0063     model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1));
-0064 end
-0065 if isfield(model,'inchis')
-0066     model.inchis(metIdx) = model.inchis(repIdx(1));
+0034 if nargin<4 || isempty(verbose)
+0035     verbose=false;
+0036 end
+0037 if nargin<5
+0038     identifiers = false;
+0039 end
+0040 
+0041 % Find occurence of replacement metabolites. Annotation will be taken from
+0042 % first metabolite found.
+0043 if identifiers
+0044     repIdx = find(strcmp(replacement,model.mets));
+0045 else
+0046     repIdx = find(strcmp(replacement,model.metNames));
+0047 end
+0048 if isempty(repIdx)
+0049     error('The replacement metabolite cannot be found in the model.');
+0050 end
+0051 
+0052 % Change name and information from metabolite to replacement metabolite
+0053 if identifiers
+0054     metIdx = find(strcmp(metabolite,model.mets));
+0055 else
+0056     metIdx = find(strcmp(metabolite,model.metNames));
+0057 end
+0058 if isempty(metIdx)
+0059     error('The to-be-replaced metabolite cannot be found in the model.');
+0060 end
+0061 
+0062 rxnsWithMet = find(model.S(metIdx,:));
+0063 if verbose==true
+0064     fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n')
+0065     fprintf(strjoin(model.rxns(rxnsWithMet),'\n'))
+0066     fprintf('\n')
 0067 end
-0068 if isfield(model,'metSmiles')
-0069     model.metSmiles(metIdx) = model.metSmiles(repIdx(1));
-0070 end
-0071 % Run through replacement metabolites and their compartments. If any of the
-0072 % to-be-replaced metabolites is already present (checked by
-0073 % metaboliteName[compartment], then the replacement metabolite is kept and
-0074 % the to-be-replace metabolite ID deleted.
-0075 
-0076 % Build list of metaboliteName[compartment]
-0077 metCompsN =cellstr(num2str(model.metComps));
-0078 map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
-0079 metCompsN = map.values(metCompsN);
-0080 metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');
-0081 
-0082 idxDelete=[];
-0083 for i = 1:length(repIdx)
-0084     metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
-0085     if length(metCompsNidx)>1
-0086         for j = 2:length(metCompsNidx)
-0087             model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
-0088             idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
-0089         end
-0090     end
-0091 end
-0092 
-0093 if ~isempty(idxDelete)
-0094     model.S(idxDelete,:) =[];
-0095     model.mets(idxDelete) = [];
-0096     model.metNames(idxDelete) = [];
-0097     model.metComps(idxDelete) = [];
-0098     model.b(idxDelete) = [];
-0099     if isfield(model,'metFormulas')
-0100         model.metFormulas(idxDelete) = [];
-0101     end
-0102     if isfield(model,'unconstrained')
-0103         model.unconstrained(idxDelete) = [];
-0104     end
-0105     if isfield(model,'metMiriams')
-0106         model.metMiriams(idxDelete) = [];
-0107     end
-0108     if isfield(model,'metCharges')
-0109         model.metCharges(idxDelete) = [];
-0110     end
-0111     if isfield(model,'metDeltaG')
-0112         model.metDeltaG(idxDelete) = [];
-0113     end
-0114     if isfield(model,'inchis')
-0115         model.inchis(idxDelete) = [];
-0116     end
-0117     if isfield(model,'metSmiles')
-0118         model.metSmiles(idxDelete) = [];
-0119     end
-0120     if isfield(model,'metFrom')
-0121         model.metFrom(idxDelete) = [];
-0122     end
-0123 end
-0124 
-0125 % This could now have created duplicate reactions. Contract model.
-0126 model=contractModel(model,[],repIdx);
-0127 end
    
+0068 
+0069 model.metNames(metIdx) = model.metNames(repIdx(1));
+0070 if isfield(model,'metFormulas')
+0071     model.metFormulas(metIdx) = model.metFormulas(repIdx(1));
+0072 end
+0073 if isfield(model,'metMiriams')
+0074     model.metMiriams(metIdx) = model.metMiriams(repIdx(1));
+0075 end
+0076 if isfield(model,'metCharges')
+0077     model.metCharges(metIdx) = model.metCharges(repIdx(1));
+0078 end
+0079 if isfield(model,'metDeltaG')
+0080     model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1));
+0081 end
+0082 if isfield(model,'inchis')
+0083     model.inchis(metIdx) = model.inchis(repIdx(1));
+0084 end
+0085 if isfield(model,'metSmiles')
+0086     model.metSmiles(metIdx) = model.metSmiles(repIdx(1));
+0087 end
+0088 
+0089 idxDelete=[];
+0090 if identifiers
+0091     originalStoch = model.S(metIdx,rxnsWithMet);
+0092     model.S(repIdx,rxnsWithMet) = originalStoch;
+0093     model.S(metIdx,rxnsWithMet) = 0;
+0094     idxDelete = metIdx;
+0095 else
+0096     % Run through replacement metabolites and their compartments. If any of the
+0097     % to-be-replaced metabolites is already present (checked by
+0098     % metaboliteName[compartment], then the replacement metabolite is kept and
+0099     % the to-be-replace metabolite ID deleted.
+0100     
+0101     % Build list of metaboliteName[compartment]
+0102     metCompsN =cellstr(num2str(model.metComps));
+0103     map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
+0104     metCompsN = map.values(metCompsN);
+0105     metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');
+0106     
+0107     for i = 1:length(repIdx)
+0108         metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
+0109         if length(metCompsNidx)>1
+0110             for j = 2:length(metCompsNidx)
+0111                 model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
+0112                 idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
+0113             end
+0114         end
+0115     end
+0116 end
+0117 
+0118 if ~isempty(idxDelete)
+0119     model.S(idxDelete,:) =[];
+0120     model.mets(idxDelete) = [];
+0121     model.metNames(idxDelete) = [];
+0122     model.metComps(idxDelete) = [];
+0123     model.b(idxDelete) = [];
+0124     if isfield(model,'metFormulas')
+0125         model.metFormulas(idxDelete) = [];
+0126     end
+0127     if isfield(model,'unconstrained')
+0128         model.unconstrained(idxDelete) = [];
+0129     end
+0130     if isfield(model,'metMiriams')
+0131         model.metMiriams(idxDelete) = [];
+0132     end
+0133     if isfield(model,'metCharges')
+0134         model.metCharges(idxDelete) = [];
+0135     end
+0136     if isfield(model,'metDeltaG')
+0137         model.metDeltaG(idxDelete) = [];
+0138     end
+0139     if isfield(model,'inchis')
+0140         model.inchis(idxDelete) = [];
+0141     end
+0142     if isfield(model,'metSmiles')
+0143         model.metSmiles(idxDelete) = [];
+0144     end
+0145     if isfield(model,'metFrom')
+0146         model.metFrom(idxDelete) = [];
+0147     end
+0148 end
+0149 
+0150 % This could now have created duplicate reactions. Contract model.
+0151 model=contractModel(model,[],repIdx);
+0152 end
    
 
    Generated by m2html © 2005
    
 
 
\ No newline at end of file

From 015bc18cf5e632bef94f1226a08a1d05cf7070c1 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven 
Date: Tue, 29 Oct 2024 14:27:01 +0100
Subject: [PATCH 8/9] refactor: use getGenesFromGrRules

---
 core/changeGrRules.m              |  19 +--
 doc/core/changeGrRules.html       |  49 +++---
 doc/core/getGenesFromGrRules.html |   2 +-
 doc/io/importModel.html           | 257 ++++++++++++++----------------
 io/importModel.m                  |  27 +---
 5 files changed, 149 insertions(+), 205 deletions(-)

diff --git a/core/changeGrRules.m b/core/changeGrRules.m
index 2ab838d3..e8418c3e 100755
--- a/core/changeGrRules.m
+++ b/core/changeGrRules.m
@@ -30,18 +30,15 @@
     error('Number of rxns and grRules should be identical')
 end
 
-for i=1:length(rxns)
-    % Add genes to model
-    geneList=transpose(cell(unique(regexp(grRules{i},'[)(]*|( and )*|( or )*','split')))); % Extract individual, unique genes from the geneAssoc provided
-    geneList=geneList(~cellfun(@isempty, geneList));
-    genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes.
-    if ~isempty(genesToAdd.genes)
-        model=addGenesRaven(model,genesToAdd); % Add genes
-    end
-    
-    % Find reaction and gene indices
-    idx=getIndexes(model,rxns,'rxns');
+% Add genes to model
+geneList = getGenesFromGrRules(grRules);
+genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes.
+if ~isempty(genesToAdd.genes)
+    model=addGenesRaven(model,genesToAdd); % Add genes
 end
+    
+% Find reaction and gene indices
+idx=getIndexes(model,rxns,'rxns');
 
 % Change gene associations
 if replace==true % Replace old gene associations
diff --git a/doc/core/changeGrRules.html b/doc/core/changeGrRules.html
index 95f7d385..82d3ea2c 100644
--- a/doc/core/changeGrRules.html
+++ b/doc/core/changeGrRules.html
@@ -48,7 +48,7 @@ 
    DESCRIPTION ^CROSS-REFERENCE INFORMATION ^
    
 This function calls:
 
+
  • addGenesRaven	addGenesRaven
  • convertCharArray	convertCharArray
  • getGenesFromGrRules	getGenesFromGrRules  Extract gene list and rxnGeneMat from grRules array.
  • getIndexes	getIndexes
  • standardizeGrRules	standardizeGrRules
    • 
 This function is called by:
 
      
 
    
@@ -89,31 +89,28 @@ 
    SOURCE CODE ^'Number of rxns and grRules should be identical')
 0031 end
 0032 
-0033 for i=1:length(rxns)
-0034     % Add genes to model
-0035     geneList=transpose(cell(unique(regexp(grRules{i},'[)(]*|( and )*|( or )*','split')))); % Extract individual, unique genes from the geneAssoc provided
-0036     geneList=geneList(~cellfun(@isempty, geneList));
-0037     genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes.
-0038     if ~isempty(genesToAdd.genes)
-0039         model=addGenesRaven(model,genesToAdd); % Add genes
-0040     end
-0041     
-0042     % Find reaction and gene indices
-0043     idx=getIndexes(model,rxns,'rxns');
-0044 end
-0045 
-0046 % Change gene associations
-0047 if replace==true % Replace old gene associations
-0048     model.grRules(idx)=grRules;
-0049 else % Add gene associations, add new gene rules after 'OR'.
-0050     model.grRules(idx)=strcat('(',model.grRules(idx),') or (',grRules,')');
-0051 end
-0052 
-0053 %Fix grRules and reconstruct rxnGeneMat
-0054 [grRules,rxnGeneMat] = standardizeGrRules(model,true);
-0055 model.grRules = grRules;
-0056 model.rxnGeneMat = rxnGeneMat;
-0057 end
    +0033 % Add genes to model +0034 geneList = getGenesFromGrRules(grRules); +0035 genesToAdd.genes=setdiff(geneList,model.genes); % Only keep the genes that are not yet part of the model.genes. +0036 if ~isempty(genesToAdd.genes) +0037 model=addGenesRaven(model,genesToAdd); % Add genes +0038 end +0039 +0040 % Find reaction and gene indices +0041 idx=getIndexes(model,rxns,'rxns'); +0042 +0043 % Change gene associations +0044 if replace==true % Replace old gene associations +0045 model.grRules(idx)=grRules; +0046 else % Add gene associations, add new gene rules after 'OR'. +0047 model.grRules(idx)=strcat('(',model.grRules(idx),') or (',grRules,')'); +0048 end +0049 +0050 %Fix grRules and reconstruct rxnGeneMat +0051 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0052 model.grRules = grRules; +0053 model.rxnGeneMat = rxnGeneMat; +0054 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/core/getGenesFromGrRules.html b/doc/core/getGenesFromGrRules.html index 15f2a098..011069b8 100644 --- a/doc/core/getGenesFromGrRules.html +++ b/doc/core/getGenesFromGrRules.html @@ -58,7 +58,7 @@

    CROSS-REFERENCE INFORMATION ^ </ul> This function is called by: <ul style= -
  • checkModelStruct checkModelStruct
    • +
  • changeGrRules changeGrRules
  • checkModelStruct checkModelStruct
    • diff --git a/doc/io/importModel.html b/doc/io/importModel.html index ed93b7fc..05f75018 100644 --- a/doc/io/importModel.html +++ b/doc/io/importModel.html @@ -105,7 +105,7 @@

    CROSS-REFERENCE INFORMATION ^ <h2><a name=SUBFUNCTIONS ^

    +
  • function fieldContent=parseNote(searchString,fieldName)
  • function fieldContent=parseAnnotation(searchString,startString,midString,fieldName)
  • function miriamStruct=parseMiriam(searchString)
  • function miriam = addSBOtoMiriam(miriam,sboTerm)

    • SOURCE CODE ^

    0001 function model=importModel(fileName,removeExcMets,removePrefix,supressWarnings)
@@ -284,13 +284,13 @@ 
    SOURCE CODE ^end
 0175 
 0176     if isfield(modelSBML.compartment(i),'annotation')
-0177         compartmentMiriams{i}=parseMiriam(modelSBML.compartment(i).annotation);
+0177         compartmentMiriams{i}=parseMiriam(modelSBML.compartment(i).annotation);
 0178     else
 0179         compartmentMiriams{i}=[];
 0180     end
 0181 
 0182     if isfield(modelSBML.compartment(i),'sboTerm') && ~(modelSBML.compartment(i).sboTerm==-1)
-0183         compartmentMiriams{i} = addSBOtoMiriam(compartmentMiriams{i},modelSBML.compartment(i).sboTerm);
+0183         compartmentMiriams{i} = addSBOtoMiriam(compartmentMiriams{i},modelSBML.compartment(i).sboTerm);
 0184     end
 0185 end
 0186 
@@ -340,7 +340,7 @@ 
    SOURCE CODE ^%Get Miriam structure
 0231         if isfield(modelSBML.species(i),'annotation')
 0232             %Get Miriam info
-0233             geneMiriam=parseMiriam(modelSBML.species(i).annotation);
+0233             geneMiriam=parseMiriam(modelSBML.species(i).annotation);
 0234             geneMiriams{numel(geneMiriams)+1,1}=geneMiriam;
 0235         else
 0236             geneMiriams{numel(geneMiriams)+1,1}=[];
@@ -352,7 +352,7 @@ 
    SOURCE CODE ^%COBRA Toolbox format. For now only the SHORT NAME is loaded
 0243         %and no mapping takes place
 0244         if isfield(modelSBML.species(i),'notes')
-0245             geneShortNames{numel(geneShortNames)+1,1}=parseNote(modelSBML.species(i).notes,'SHORT NAME');
+0245             geneShortNames{numel(geneShortNames)+1,1}=parseNote(modelSBML.species(i).notes,'SHORT NAME');
 0246         else
 0247             geneShortNames{numel(geneShortNames)+1,1}='';
 0248         end
@@ -425,12 +425,12 @@ 
    SOURCE CODE ^end
 0316 
 0317             %Get Miriam info
-0318             metMiriam=parseMiriam(modelSBML.species(i).annotation);
+0318             metMiriam=parseMiriam(modelSBML.species(i).annotation);
 0319             metaboliteMiriams{numel(metaboliteMiriams)+1,1}=metMiriam;
 0320         else
 0321             metaboliteInChI{numel(metaboliteInChI)+1,1}='';
 0322             if isfield(modelSBML.species(i),'notes')
-0323                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
+0323                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
 0324             else
 0325                 metaboliteFormula{numel(metaboliteFormula)+1,1}='';
 0326             end
@@ -438,7 +438,7 @@ 
    SOURCE CODE ^end
 0329         if ~isempty(modelSBML.species(i).notes)
 0330             if ~isfield(modelSBML.species(i),'annotation')
-0331                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
+0331                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
 0332             end
 0333         elseif ~isfield(modelSBML.species(i),'annotation')
 0334             metaboliteFormula{numel(metaboliteFormula)+1,1}='';
@@ -461,7 +461,7 @@ 
    SOURCE CODE ^else
 0352                     if isfield(modelSBML.species(i),'notes')
 0353                         if strfind(modelSBML.species(i).notes,'CHARGE')
-0354                             metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
+0354                             metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
 0355                         else
 0356                             metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
 0357                         end
@@ -471,7 +471,7 @@ 
    SOURCE CODE ^end
 0362             elseif isfield(modelSBML.species(i),'notes')
 0363                 if strfind(modelSBML.species(i).notes,'CHARGE')
-0364                     metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
+0364                     metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
 0365                 else
 0366                     metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
 0367                 end
@@ -491,12 +491,12 @@ 
    SOURCE CODE ^%Add SBO terms to gene and metabolite miriam fields
 0382 if numel(unique(geneSBOs)) > 1  % don't add if they're all identical
 0383     for i = 1:numel(geneNames)
-0384         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},geneSBOs(i));
+0384         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},geneSBOs(i));
 0385     end
 0386 end
 0387 if numel(unique(metSBOs)) > 1
 0388     for i = 1:numel(metaboliteNames)
-0389         metaboliteMiriams{i} = addSBOtoMiriam(metaboliteMiriams{i},metSBOs(i));
+0389         metaboliteMiriams{i} = addSBOtoMiriam(metaboliteMiriams{i},metSBOs(i));
 0390     end
 0391 end
 0392 
@@ -603,9 +603,9 @@ 
    SOURCE CODE ^%it will be executed, if 'GENE_ASSOCIATION' is found in
 0494         %modelSBML.reaction(i).notes
 0495         if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION')
-0496             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION');
+0496             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION');
 0497         elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION')
-0498             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION');
+0498             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION');
 0499         else
 0500             geneAssociation='';
 0501         end
@@ -677,7 +677,7 @@ 
    SOURCE CODE ^'';
 0568         end
 0569     elseif isfield(modelSBML.reaction(i),'notes')
-0570         rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT');
+0570         rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT');
 0571     end
 0572     if ~isempty(rxnComp)
 0573         %Find it in the compartment list
@@ -686,42 +686,42 @@ 
    SOURCE CODE ^end
 0577 
 0578 
-0579     miriamStruct=parseMiriam(modelSBML.reaction(i).annotation);
+0579     miriamStruct=parseMiriam(modelSBML.reaction(i).annotation);
 0580     rxnMiriams{counter}=miriamStruct;
 0581     if isfield(modelSBML.reaction(i),'notes')
-0582         subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM'));
+0582         subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM'));
 0583         subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = [];
 0584         if strfind(modelSBML.reaction(i).notes,'Confidence Level')
-0585             confScore = parseNote(modelSBML.reaction(i).notes,'Confidence Level');
+0585             confScore = parseNote(modelSBML.reaction(i).notes,'Confidence Level');
 0586             if isempty(confScore)
 0587                 confScore = 0;
 0588             end
 0589             rxnconfidencescores(counter)=str2double(confScore);
 0590         end
-0591         rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS');
-0592         rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES');
+0591         rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS');
+0592         rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES');
 0593     end
 0594 
 0595     %Get SBO terms
 0596     if isfield(modelSBML.reaction(i),'sboTerm') && ~(modelSBML.reaction(i).sboTerm==-1)
-0597         rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm);
+0597         rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm);
 0598     end
 0599 
 0600     %Get ec-codes
 0601     eccode='';
 0602     if ~isempty(modelSBML.reaction(i).annotation)
 0603         if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code')
-0604             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code');
+0604             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code');
 0605         elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code')
-0606             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code');
+0606             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code');
 0607         elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code')
-0608             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code');
+0608             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code');
 0609         end
 0610     elseif isfield(modelSBML.reaction(i),'notes')
 0611         if strfind(modelSBML.reaction(i).notes,'EC Number')
-0612             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')];
+0612             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')];
 0613         elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS')
-0614             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')];
+0614             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')];
 0615         end
 0616     end
 0617     eccodes{counter}=eccode;
@@ -969,15 +969,15 @@ 
    SOURCE CODE ^'sboTerm');
 0860                 end
 0861                 for i = 1:numel(genes)
-0862                     geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation);
+0862                     geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation);
 0863                     if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1)
-0864                         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm);
+0864                         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm);
 0865                     end
 0866                 end
 0867             end
 0868             proteins={modelSBML.fbc_geneProduct.fbc_name};
 0869         else
-0870             genes=getGeneList(grRules);
+0870             genes=getGenesFromGrRules(grRules);
 0871         end
 0872         model.genes=genes;
 0873         model.grRules=grRules;
@@ -1169,122 +1169,97 @@ 
    SOURCE CODE ^end
 1060 end
 1061 
-1062 function matchGenes=getGeneList(grRules)
-1063 %Constructs the list of unique genes from grRules
-1064 
-1065 %Assumes that everything that isn't a paranthesis, " AND " or " or " is a
-1066 %gene name
-1067 genes=strrep(grRules,'(','');
-1068 genes=strrep(genes,')','');
-1069 genes=strrep(genes,' or ',' ');
-1070 genes=strrep(genes,' and ',' ');
-1071 genes=strrep(genes,' OR ',' ');
-1072 genes=strrep(genes,' AND ',' ');
-1073 genes=regexp(genes,' ','split');
-1074 
-1075 allNames={};
-1076 for i=1:numel(genes)
-1077     allNames=[allNames genes{i}];
+1062 function fieldContent=parseNote(searchString,fieldName)
+1063 %The function obtains the particular information from 'notes' field, using
+1064 %fieldName as the dummy string
+1065 
+1066 fieldContent='';
+1067 
+1068 if strfind(searchString,fieldName)
+1069     [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match');
+1070     targetString=regexprep(targetString,'<p>|</p>','');
+1071     targetString=regexprep(targetString,[fieldName, ':'],'');
+1072     for i=1:numel(targetString)
+1073         fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
+1074     end
+1075     fieldContent=regexprep(fieldContent,'^;|;$','');
+1076 else
+1077     fieldContent='';
 1078 end
-1079 matchGenes=unique(allNames)';
+1079 end
 1080 
-1081 %Remove the empty element if present
-1082 if isempty(matchGenes{1})
-1083     matchGenes(1)=[];
-1084 end
-1085 end
-1086 
-1087 function fieldContent=parseNote(searchString,fieldName)
-1088 %The function obtains the particular information from 'notes' field, using
-1089 %fieldName as the dummy string
-1090 
-1091 fieldContent='';
+1081 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName)
+1082 
+1083 fieldContent='';
+1084 
+1085 %Removing whitespace characters from the ending strings, which may occur in
+1086 %several cases
+1087 searchString=regexprep(searchString,'" />','"/>');
+1088 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match');
+1089 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
+1090 targetString=regexprep(targetString,startString,'');
+1091 targetString=regexprep(targetString,[fieldName midString],'');
 1092 
-1093 if strfind(searchString,fieldName)
-1094     [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match');
-1095     targetString=regexprep(targetString,'<p>|</p>','');
-1096     targetString=regexprep(targetString,[fieldName, ':'],'');
-1097     for i=1:numel(targetString)
-1098         fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
-1099     end
-1100     fieldContent=regexprep(fieldContent,'^;|;$','');
-1101 else
-1102     fieldContent='';
-1103 end
-1104 end
-1105 
-1106 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName)
-1107 
-1108 fieldContent='';
-1109 
-1110 %Removing whitespace characters from the ending strings, which may occur in
-1111 %several cases
-1112 searchString=regexprep(searchString,'" />','"/>');
-1113 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match');
-1114 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
-1115 targetString=regexprep(targetString,startString,'');
-1116 targetString=regexprep(targetString,[fieldName midString],'');
-1117 
-1118 for i=1:numel(targetString)
-1119     fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
-1120 end
+1093 for i=1:numel(targetString)
+1094     fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
+1095 end
+1096 
+1097 fieldContent=regexprep(fieldContent,'^;|;$','');
+1098 end
+1099 
+1100 function miriamStruct=parseMiriam(searchString)
+1101 %Generates miriam structure from annotation field
+1102 
+1103 %Finding whether miriams are written in the old or the new way
+1104 if strfind(searchString,'urn:miriam:')
+1105     startString='urn:miriam:';
+1106 elseif strfind(searchString,'http://identifiers.org/')
+1107     startString='http://identifiers.org/';
+1108 elseif strfind(searchString,'https://identifiers.org/')
+1109     startString='https://identifiers.org/';
+1110 else
+1111     miriamStruct=[];
+1112     return;
+1113 end
+1114 
+1115 miriamStruct=[];
+1116 
+1117 searchString=regexprep(searchString,'" />','"/>');
+1118 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match');
+1119 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
+1120 targetString=regexprep(targetString,startString,'');
 1121 
-1122 fieldContent=regexprep(fieldContent,'^;|;$','');
-1123 end
-1124 
-1125 function miriamStruct=parseMiriam(searchString)
-1126 %Generates miriam structure from annotation field
-1127 
-1128 %Finding whether miriams are written in the old or the new way
-1129 if strfind(searchString,'urn:miriam:')
-1130     startString='urn:miriam:';
-1131 elseif strfind(searchString,'http://identifiers.org/')
-1132     startString='http://identifiers.org/';
-1133 elseif strfind(searchString,'https://identifiers.org/')
-1134     startString='https://identifiers.org/';
-1135 else
-1136     miriamStruct=[];
-1137     return;
-1138 end
-1139 
-1140 miriamStruct=[];
-1141 
-1142 searchString=regexprep(searchString,'" />','"/>');
-1143 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match');
-1144 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
-1145 targetString=regexprep(targetString,startString,'');
-1146 
-1147 fwdslash  = contains(targetString,'/');
-1148 midString = cell(numel(targetString),1);
-1149 midString(fwdslash) = {'/'};
-1150 midString(~fwdslash) = {':'};
-1151 
-1152 counter=0;
-1153 for i=1:numel(targetString)
-1154     if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once'))
-1155         counter=counter+1;
-1156         miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once');
-1157         miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once');
-1158         miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.','');
-1159     end
-1160 end
-1161 end
-1162 
-1163 function miriam = addSBOtoMiriam(miriam,sboTerm)
-1164 %Appends SBO term to miriam structure
-1165 
-1166 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]};  % convert to proper format
-1167 if isempty(miriam)
-1168     miriam.name = {'sbo'};
-1169     miriam.value = sboTerm;
-1170 elseif any(strcmp('sbo',miriam.name))
-1171     currSbo = strcmp('sbo',miriam.name);
-1172     miriam.value(currSbo) = sboTerm;
-1173 else
-1174     miriam.name(end+1) = {'sbo'};
-1175     miriam.value(end+1) = sboTerm;
-1176 end
-1177 end
    +1122 fwdslash = contains(targetString,'/'); +1123 midString = cell(numel(targetString),1); +1124 midString(fwdslash) = {'/'}; +1125 midString(~fwdslash) = {':'}; +1126 +1127 counter=0; +1128 for i=1:numel(targetString) +1129 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) +1130 counter=counter+1; +1131 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); +1132 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); +1133 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); +1134 end +1135 end +1136 end +1137 +1138 function miriam = addSBOtoMiriam(miriam,sboTerm) +1139 %Appends SBO term to miriam structure +1140 +1141 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format +1142 if isempty(miriam) +1143 miriam.name = {'sbo'}; +1144 miriam.value = sboTerm; +1145 elseif any(strcmp('sbo',miriam.name)) +1146 currSbo = strcmp('sbo',miriam.name); +1147 miriam.value(currSbo) = sboTerm; +1148 else +1149 miriam.name(end+1) = {'sbo'}; +1150 miriam.value(end+1) = sboTerm; +1151 end +1152 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/io/importModel.m b/io/importModel.m index 3032ba53..e8a412d4 100755 --- a/io/importModel.m +++ b/io/importModel.m @@ -867,7 +867,7 @@ end proteins={modelSBML.fbc_geneProduct.fbc_name}; else - genes=getGeneList(grRules); + genes=getGenesFromGrRules(grRules); end model.genes=genes; model.grRules=grRules; @@ -1059,31 +1059,6 @@ end end -function matchGenes=getGeneList(grRules) -%Constructs the list of unique genes from grRules - -%Assumes that everything that isn't a paranthesis, " AND " or " or " is a -%gene name -genes=strrep(grRules,'(',''); -genes=strrep(genes,')',''); -genes=strrep(genes,' or ',' '); -genes=strrep(genes,' and ',' '); -genes=strrep(genes,' OR ',' '); -genes=strrep(genes,' AND ',' '); -genes=regexp(genes,' ','split'); - -allNames={}; -for i=1:numel(genes) - allNames=[allNames genes{i}]; -end -matchGenes=unique(allNames)'; - -%Remove the empty element if present -if isempty(matchGenes{1}) - matchGenes(1)=[]; -end -end - function fieldContent=parseNote(searchString,fieldName) %The function obtains the particular information from 'notes' field, using %fieldName as the dummy string From d90bce85f72423f56f1b65a5c0d66bcbcb8ae941 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Thu, 14 Nov 2024 14:01:19 +0100 Subject: [PATCH 9/9] doc: checkRxn clarify meaning of output --- core/checkRxn.m | 8 ++++---- doc/core/checkRxn.html | 16 ++++++++-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/core/checkRxn.m b/core/checkRxn.m index cf29bf94..1334514c 100755 --- a/core/checkRxn.m +++ b/core/checkRxn.m @@ -14,11 +14,11 @@ % % report % reactants array with reactant indexes -% canMake boolean array, true if the corresponding reactant can be -% synthesized +% canMake boolean array, true if the corresponding reactant can +% be synthesized by the rest of the metabolic network % products array with product indexes -% canConsume boolean array, true if the corresponding reactant can -% be consumed +% canConsume boolean array, true if the corresponding product can +% be consumed by the rest of the metabolic network % % Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport) diff --git a/doc/core/checkRxn.html b/doc/core/checkRxn.html index 2c31271b..4169caa7 100644 --- a/doc/core/checkRxn.html +++ b/doc/core/checkRxn.html @@ -42,11 +42,11 @@

    DESCRIPTION ^SOURCE CODE ^% 0015 % report 0016 % reactants array with reactant indexes -0017 % canMake boolean array, true if the corresponding reactant can be -0018 % synthesized +0017 % canMake boolean array, true if the corresponding reactant can +0018 % be synthesized by the rest of the metabolic network 0019 % products array with product indexes -0020 % canConsume boolean array, true if the corresponding reactant can -0021 % be consumed +0020 % canConsume boolean array, true if the corresponding product can +0021 % be consumed by the rest of the metabolic network 0022 % 0023 % Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport) 0024