From 5b6e3c64120b1260fdc575acda4a8c7556c6a350 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 7 Jul 2026 15:48:01 +0200 Subject: [PATCH 1/2] fix: gene field alignment in removeLowScoreGenes removeLowScoreGenes regenerated model.genes from getGenesFromGrRules, which returns a sorted list, but trimmed the gene-associated fields (geneShortNames, proteins, geneMiriams, geneFrom, geneComps) using a mask in the original gene order. When model.genes was not already sorted, this left every annotation field shifted relative to model.genes, corrupting the gene ID to gene symbol mapping in ftINIT-reconstructed models (removeGenes=true). Retain the remaining genes in their original order and reorder the rxnGeneMat columns to match, so the trimmed annotation fields stay aligned. Fixes #669 --- INIT/removeLowScoreGenes.m | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/INIT/removeLowScoreGenes.m b/INIT/removeLowScoreGenes.m index 76a5bdeb..a98a65cf 100644 --- a/INIT/removeLowScoreGenes.m +++ b/INIT/removeLowScoreGenes.m @@ -109,13 +109,20 @@ % regenerate "genes" and "rxnGeneMat" model fields [genes,rxnGeneMat] = getGenesFromGrRules(newModel.grRules); -newModel.genes = genes; -newModel.rxnGeneMat = rxnGeneMat; -% update other gene-related fields -remInd = ~ismember(model.genes,newModel.genes); +% determine which of the original genes were removed +remInd = ~ismember(model.genes,genes); remGenes = model.genes(remInd); +% Keep the retained genes in their original order rather than the sorted +% order returned by getGenesFromGrRules. Gene removal never introduces new +% genes, so the remaining genes are model.genes(~remInd). Preserving this +% order ensures the gene-associated fields trimmed below (geneShortNames, +% proteins, etc.), which are indexed by remInd, stay aligned with genes. +newModel.genes = model.genes(~remInd); +[~,reorderInd] = ismember(newModel.genes,genes); +newModel.rxnGeneMat = rxnGeneMat(:,reorderInd); + if isfield(newModel,'geneShortNames') newModel.geneShortNames(remInd) = []; end From aed0be308c7c9b2eb3848b5688897f79744e2016 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 7 Jul 2026 15:53:11 +0200 Subject: [PATCH 2/2] chore: updateDocumentation --- doc/INIT/INITStepDesc.html | 100 ++-- doc/INIT/ftINIT.html | 407 ++++++++-------- doc/INIT/ftINITFillGapsMILP.html | 49 +- doc/INIT/removeLowScoreGenes.html | 297 ++++++------ doc/io/exportModel.html | 451 +++++++++--------- doc/testing/unit_tests/importExportTests.html | 56 ++- 6 files changed, 705 insertions(+), 655 deletions(-) diff --git a/doc/INIT/INITStepDesc.html b/doc/INIT/INITStepDesc.html index 02850578..113533f7 100644 --- a/doc/INIT/INITStepDesc.html +++ b/doc/INIT/INITStepDesc.html @@ -70,55 +70,57 @@

SOURCE CODE ^% .mets Names of metabolites to remove 0025 % .compsToKeep Compartments for which metabolites should be kept. 0026 MILPParams %Cell array of MILPparams - dictates how many iterations that will be run in this step. -0027 %Typically, MIPGap and TimeLimit is specified -0028 AbsMIPGaps %If the objective is close to zero, a percentage of that is very small. -0029 %Therefore, also set an absolut value for this (typically 10 or 20). -0030 %For practical reasons, the first number is not used -0031 end -0032 methods -0033 function obj = INITStepDesc(posRevOff_, AllowMetSecr_, howToUsePrevResults_, rxnsToIgnoreMask_, metsToIgnore_, MILPParams_, absMIPGaps_) -0034 if nargin > 0 -0035 obj.PosRevOff = posRevOff_; -0036 else -0037 obj.PosRevOff = false; -0038 end -0039 if nargin > 1 -0040 obj.AllowMetSecr = AllowMetSecr_; -0041 else -0042 obj.AllowMetSecr = false; -0043 end -0044 if nargin > 2 -0045 obj.HowToUsePrevResults = howToUsePrevResults_; -0046 else -0047 obj.HowToUsePrevResults = 'essential'; -0048 end -0049 if nargin > 3 -0050 obj.RxnsToIgnoreMask = rxnsToIgnoreMask_; -0051 else -0052 obj.RxnsToIgnoreMask = [1;0;0;0;0;0;0;0]; -0053 end -0054 if nargin > 4 -0055 obj.MetsToIgnore = metsToIgnore_; -0056 else -0057 obj.MetsToIgnore = [1;0;0;0;0;0;0;0]; -0058 end -0059 if nargin > 5 -0060 obj.MILPParams = MILPParams_; -0061 else -0062 params = struct(); -0063 params.TimeLimit = 5000; -0064 params.MIPGap = 0.0004; -0065 obj.MILPParams = {params}; -0066 end -0067 -0068 if nargin > 6 -0069 obj.AbsMIPGaps = absMIPGaps_; -0070 else -0071 obj.AbsMIPGaps = 10; -0072 end -0073 end -0074 end -0075 end +0027 %Typically, MIPGap and TimeLimit is specified. ftINIT additionally +0028 %defaults Threads to 1 (single-threaded Gurobi) for determinism; +0029 %set Threads here to override (0 = use all cores). +0030 AbsMIPGaps %If the objective is close to zero, a percentage of that is very small. +0031 %Therefore, also set an absolut value for this (typically 10 or 20). +0032 %For practical reasons, the first number is not used +0033 end +0034 methods +0035 function obj = INITStepDesc(posRevOff_, AllowMetSecr_, howToUsePrevResults_, rxnsToIgnoreMask_, metsToIgnore_, MILPParams_, absMIPGaps_) +0036 if nargin > 0 +0037 obj.PosRevOff = posRevOff_; +0038 else +0039 obj.PosRevOff = false; +0040 end +0041 if nargin > 1 +0042 obj.AllowMetSecr = AllowMetSecr_; +0043 else +0044 obj.AllowMetSecr = false; +0045 end +0046 if nargin > 2 +0047 obj.HowToUsePrevResults = howToUsePrevResults_; +0048 else +0049 obj.HowToUsePrevResults = 'essential'; +0050 end +0051 if nargin > 3 +0052 obj.RxnsToIgnoreMask = rxnsToIgnoreMask_; +0053 else +0054 obj.RxnsToIgnoreMask = [1;0;0;0;0;0;0;0]; +0055 end +0056 if nargin > 4 +0057 obj.MetsToIgnore = metsToIgnore_; +0058 else +0059 obj.MetsToIgnore = [1;0;0;0;0;0;0;0]; +0060 end +0061 if nargin > 5 +0062 obj.MILPParams = MILPParams_; +0063 else +0064 params = struct(); +0065 params.TimeLimit = 5000; +0066 params.MIPGap = 0.0004; +0067 obj.MILPParams = {params}; +0068 end +0069 +0070 if nargin > 6 +0071 obj.AbsMIPGaps = absMIPGaps_; +0072 else +0073 obj.AbsMIPGaps = 10; +0074 end +0075 end +0076 end +0077 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/INIT/ftINIT.html b/doc/INIT/ftINIT.html index 3c05d6e0..6133e894 100644 --- a/doc/INIT/ftINIT.html +++ b/doc/INIT/ftINIT.html @@ -368,206 +368,213 @@

SOURCE CODE ^if ~isfield(params, 'TimeLimit') 0234 params.TimeLimit = 5000; 0235 end -0236 -0237 if ~first -0238 %There is sometimes a problem with that the objective function becomes close to zero, -0239 %which leads to that a small percentage of that (which is the MIPGap sent in) is very small -0240 %and the MILP hence takes a lot of time to finish. We also therefore use an absolute MIP gap, -0241 %converted to a percentage using the last value of the objective function. -0242 params.MIPGap = min(max(params.MIPGap, stp.AbsMIPGaps{rn}/abs(lastObjVal)),1); -0243 params.seed = 1234;%use another seed, may work better -0244 -0245 if mipGap <= params.MIPGap -0246 success = true; -0247 break; %we're done - this will not happen the first time -0248 else -0249 disp(['MipGap too high, trying with a different run. MipGap = ' num2str(mipGap) ' New MipGap Limit = ' num2str(params.MIPGap)]) -0250 end -0251 end -0252 -0253 first = false; -0254 -0255 %now run the MILP -0256 try -0257 %The prodweight for metabolomics is currently set to 5 - 0.5 was default in the old version, which I deemed very small? -0258 %There could be a need to specify this somewhere in the call at some point. -0259 %This value has not been evaluated, but is assumed in the test cases - if changed, update the test case -0260 startVals = []; -0261 if ~isempty(fullMipRes) -0262 startVals = fullMipRes.full; -0263 end -0264 [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose); -0265 %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed -0266 fluxes1(toRev) = -fluxes1(toRev); -0267 -0268 mipGap = fullMipRes.mipgap; -0269 lastObjVal = fullMipRes.obj; -0270 catch e -0271 mipGap = Inf; -0272 lastObjVal = Inf; %we need to set something here, Inf leads to that this doesn't come into play -0273 end -0274 -0275 success = mipGap <= params.MIPGap; -0276 end -0277 -0278 if ~success -0279 dispEM(['Failed to find good enough solution within the time frame. MIPGap: ' num2str(mipGap)]); -0280 end -0281 -0282 %save the reactions turned on and their fluxes for the next step -0283 rxnsTurnedOn = rxnsTurnedOn | rxnsTurnedOn1.'; -0284 %The fluxes are a bit tricky - what if they change direction between the steps? -0285 %The fluxes are used to determine the direction in which reactions are forced on -0286 %(to simplify the problem it is good if they are unidirectional). -0287 %We use the following strategy: -0288 %1. Use the fluxes from the most recent step. -0289 %2. If any flux is very low there (i.e. basically zero), use the flux from the previous steps -0290 %This could in theory cause problems, but seems to work well practically -0291 fluxesOld = fluxes; -0292 fluxes = fluxes1; -0293 %make sure that all reactions that are on actually has a flux - otherwise -0294 %things could go bad, since the flux will be set to essential in a random direction -0295 %This sometimes happens for rxns with negative score - let's just accept that. -0296 %if (sum(abs(fluxes1) < 10^-7 & rxnsTurnedOn)) -0297 % dispEM('There are rxns turned on without flux - this might cause problems'); -0298 %end -0299 %fluxes(abs(fluxes1) < 10^-7) = fluxesOld(abs(fluxes1) < 10^-9); -0300 end -0301 -0302 -0303 %get the essential rxns -0304 essential = ismember(prepData.minModel.rxns,prepData.essentialRxns); -0305 %So, we only add reactions where the linearly merged scores are zero for all linearly dependent reactions -0306 % (this cannot happen by chance, taken care of in the function groupRxnScores) -0307 rxnsToIgn = rxnScores == 0; -0308 deletedRxnsInINITSel = ~(rxnsTurnedOn | rxnsToIgn | essential); -0309 deletedRxnsInINIT = prepData.minModel.rxns(deletedRxnsInINITSel); -0310 -0311 %Here we need to figure out which original reactions (before the linear merge) -0312 %that were removed. These are all reactions with the same group ids as the removed reactions -0313 groupIdsRemoved = prepData.groupIds(ismember(prepData.refModel.rxns, deletedRxnsInINIT)); %can improve this slightly, use sel above -0314 groupIdsRemoved = groupIdsRemoved(groupIdsRemoved ~= 0);%zero means that the reaction was not grouped, all with zeros are not a group! -0315 rxnsToRem = union(prepData.refModel.rxns(ismember(prepData.groupIds,groupIdsRemoved)), deletedRxnsInINIT);%make a union here to include the ungrouped (unmerged) as well -0316 -0317 initModel = removeReactions(prepData.refModel,rxnsToRem,false,true); -0318 -0319 % remove metabolites separately to avoid removing those needed for tasks -0320 unusedMets = initModel.mets(all(initModel.S == 0,2)); -0321 initModel = removeMets(initModel, setdiff(unusedMets, prepData.essentialMetsForTasks)); -0322 -0323 %if printReport == true -0324 % printScores(initModel,'INIT model statistics',hpaData,transcrData,tissue,celltype); -0325 % printScores(removeReactions(cModel,setdiff(cModel.rxns,rxnsToRem),true,true),'Reactions deleted by INIT',hpaData,transcrData,tissue,celltype); -0326 %end -0327 -0328 %The full model has exchange reactions in it. ftINITFillGapsForAllTasks calls -0329 %ftINITFillGaps, which automatically removes exchange metabolites (because it -0330 %assumes that the reactions are constrained when appropriate). In this case the -0331 %uptakes/outputs are retrieved from the task sheet instead. To prevent -0332 %exchange reactions being used to fill gaps, they are deleted from the -0333 %reference model here. -0334 initModel.id = 'INITModel'; -0335 -0336 %If gaps in the model should be filled using a task list -0337 if ~isempty(prepData.taskStruct) -0338 %Remove exchange reactions and reactions already included in the INIT -0339 %model -0340 %We changed strategy and instead include all rxns except the exchange rxns in the ref model -0341 %But we do keep the exchange rxns that are essential. -0342 %Let's test to remove all, that should work -0343 -0344 %At this stage the model is fully connected and most of the genes with -0345 %good scores should have been included. The final gap-filling should -0346 %take the scores of the genes into account, so that "rather bad" -0347 %reactions are preferred to "very bad" reactions. However, reactions -0348 %with positive scores will be included even if they are not connected -0349 %in the current formulation. Therefore, such reactions will have to be -0350 %assigned a small negative score instead. -0351 exchRxns = getExchangeRxns(prepData.refModel); -0352 refModelNoExc = removeReactions(prepData.refModelWithBM,exchRxns,false,true); -0353 exchRxns = getExchangeRxns(initModel); -0354 initModelNoExc = removeReactions(closeModel(initModel),exchRxns,false,true); -0355 -0356 if useScoresForTasks == true -0357 %map the rxn scores to the model without exchange rxns -0358 [~,ia,ib] = intersect(refModelNoExc.rxns,prepData.refModel.rxns); -0359 rxnScores2nd = NaN(length(refModelNoExc.rxns),1); -0360 rxnScores2nd(ia) = origRxnScores(ib); -0361 %all(rxnScores2nd == refRxnScores);%should be the same, ok! -0362 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,min(rxnScores2nd,-0.1),prepData.taskStruct,paramsFT,verbose); -0363 else -0364 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,[],prepData.taskStruct,paramsFT,verbose); -0365 end -0366 %if printReport == true -0367 % printScores(outModel,'Functional model statistics',hpaData,transcrData,tissue,celltype); -0368 % printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,transcrData,tissue,celltype); -0369 %end -0370 -0371 addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2)); -0372 else -0373 outModel = initModel; -0374 addedRxnMat = []; -0375 addedRxnsForTasks = {}; -0376 end -0377 -0378 % The model can now perform all the tasks defined in the task list. -0379 model = outModel; -0380 -0381 -0382 % At this stage the model will contain some exchange reactions but probably -0383 % not all (and maybe zero). This can be inconvenient, so all exchange -0384 % reactions from the reference model are added, except for those which -0385 % involve metabolites that are not in the model. -0386 -0387 %Start from the original model, and just remove the reactions that are no longer there (and keep exchange rxns). The model we got out -0388 %from the problem is not complete, it doesn't have GRPs etc. -0389 %The logic below is a bit complicated. We identify the reactions that should be removed from the full model as -0390 %reactions that have been removed in the init model except the ones that were added back. In addition, we make -0391 %sure that no exchange rxns are removed - they can be removed in the init model if they were linearly merged with other -0392 %reactions that were decided to be removed from the model. We want to keep all exchange rxns to make sure the tasks can -0393 %be performed also without manipulating the b vector in the model (which is what is done in the gap-filling). -0394 exchRxns = getExchangeRxns(prepData.refModel); -0395 deletedRxnsInINIT = setdiff(prepData.refModel.rxns,union(union(initModel.rxns, addedRxnsForTasks), exchRxns)); -0396 outModel = removeReactions(prepData.refModel, deletedRxnsInINIT, true); %we skip removing the genes for now, I'm not sure it is desirable -0397 -0398 % If requested, attempt to remove negative-score genes from the model, -0399 % depending on their role (isozyme or complex subunit) in each grRule. -0400 % See the "removeLowScoreGenes" function more more details, and to adjust -0401 % any default parameters therein. -0402 if ( removeGenes ) -0403 [~, geneScores] = scoreComplexModel(outModel,hpaData,transcrData,tissue,celltype); -0404 outModel = removeLowScoreGenes(outModel,geneScores); -0405 end -0406 -0407 -0408 model = outModel; -0409 -0410 end -0411 -0412 %This is for printing a summary of a model -0413 function [rxnS, geneS] = printScores(model,name,hpaData,transcrData,tissue,celltype) -0414 [a, b] = scoreComplexModel(model,hpaData,transcrData,tissue,celltype); -0415 rxnS = mean(a); -0416 geneS = mean(b,'omitnan'); -0417 fprintf([name ':\n']); -0418 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']); -0419 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']); -0420 fprintf(['\tMean gene score: ' num2str(geneS) '\n']); -0421 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']); -0422 end -0423 -0424 function rxnsToIgnore = getRxnsFromPattern(rxnsToIgnorePattern, prepData) -0425 rxnsToIgnore = false(length(prepData.toIgnoreExch),1); -0426 if rxnsToIgnorePattern(1) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreExch; end; -0427 if rxnsToIgnorePattern(2) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreImportRxns; end; -0428 if rxnsToIgnorePattern(3) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSimpleTransp; end; -0429 if rxnsToIgnorePattern(4) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAdvTransp; end; -0430 if rxnsToIgnorePattern(5) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSpont; end; -0431 if rxnsToIgnorePattern(6) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreS; end; -0432 if rxnsToIgnorePattern(7) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreCustomRxns; end; -0433 if rxnsToIgnorePattern(8) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAllWithoutGPRs; end; -0434 end -0435 +0236 +0237 %Default to single-threaded Gurobi MILP solving. Multi-threaded Gurobi +0238 %can non-deterministically report the MILP as infeasible (issue #607). +0239 %Override by setting 'Threads' in the step's MILPParams (0 = all cores). +0240 if ~isfield(params, 'Threads') +0241 params.Threads = 1; +0242 end +0243 +0244 if ~first +0245 %There is sometimes a problem with that the objective function becomes close to zero, +0246 %which leads to that a small percentage of that (which is the MIPGap sent in) is very small +0247 %and the MILP hence takes a lot of time to finish. We also therefore use an absolute MIP gap, +0248 %converted to a percentage using the last value of the objective function. +0249 params.MIPGap = min(max(params.MIPGap, stp.AbsMIPGaps{rn}/abs(lastObjVal)),1); +0250 params.seed = 1234;%use another seed, may work better +0251 +0252 if mipGap <= params.MIPGap +0253 success = true; +0254 break; %we're done - this will not happen the first time +0255 else +0256 disp(['MipGap too high, trying with a different run. MipGap = ' num2str(mipGap) ' New MipGap Limit = ' num2str(params.MIPGap)]) +0257 end +0258 end +0259 +0260 first = false; +0261 +0262 %now run the MILP +0263 try +0264 %The prodweight for metabolomics is currently set to 5 - 0.5 was default in the old version, which I deemed very small? +0265 %There could be a need to specify this somewhere in the call at some point. +0266 %This value has not been evaluated, but is assumed in the test cases - if changed, update the test case +0267 startVals = []; +0268 if ~isempty(fullMipRes) +0269 startVals = fullMipRes.full; +0270 end +0271 [deletedRxnsInINIT1, metProduction,fullMipRes,rxnsTurnedOn1,fluxes1] = ftINITInternalAlg(mm,rxnScores,metData,essentialRxns,5,stp.AllowMetSecr,stp.PosRevOff,params, startVals, fluxes, verbose); +0272 %This is a bit tricky - since we reversed some reactions, those fluxes also need to be reversed +0273 fluxes1(toRev) = -fluxes1(toRev); +0274 +0275 mipGap = fullMipRes.mipgap; +0276 lastObjVal = fullMipRes.obj; +0277 catch e +0278 mipGap = Inf; +0279 lastObjVal = Inf; %we need to set something here, Inf leads to that this doesn't come into play +0280 end +0281 +0282 success = mipGap <= params.MIPGap; +0283 end +0284 +0285 if ~success +0286 dispEM(['Failed to find good enough solution within the time frame. MIPGap: ' num2str(mipGap)]); +0287 end +0288 +0289 %save the reactions turned on and their fluxes for the next step +0290 rxnsTurnedOn = rxnsTurnedOn | rxnsTurnedOn1.'; +0291 %The fluxes are a bit tricky - what if they change direction between the steps? +0292 %The fluxes are used to determine the direction in which reactions are forced on +0293 %(to simplify the problem it is good if they are unidirectional). +0294 %We use the following strategy: +0295 %1. Use the fluxes from the most recent step. +0296 %2. If any flux is very low there (i.e. basically zero), use the flux from the previous steps +0297 %This could in theory cause problems, but seems to work well practically +0298 fluxesOld = fluxes; +0299 fluxes = fluxes1; +0300 %make sure that all reactions that are on actually has a flux - otherwise +0301 %things could go bad, since the flux will be set to essential in a random direction +0302 %This sometimes happens for rxns with negative score - let's just accept that. +0303 %if (sum(abs(fluxes1) < 10^-7 & rxnsTurnedOn)) +0304 % dispEM('There are rxns turned on without flux - this might cause problems'); +0305 %end +0306 %fluxes(abs(fluxes1) < 10^-7) = fluxesOld(abs(fluxes1) < 10^-9); +0307 end +0308 +0309 +0310 %get the essential rxns +0311 essential = ismember(prepData.minModel.rxns,prepData.essentialRxns); +0312 %So, we only add reactions where the linearly merged scores are zero for all linearly dependent reactions +0313 % (this cannot happen by chance, taken care of in the function groupRxnScores) +0314 rxnsToIgn = rxnScores == 0; +0315 deletedRxnsInINITSel = ~(rxnsTurnedOn | rxnsToIgn | essential); +0316 deletedRxnsInINIT = prepData.minModel.rxns(deletedRxnsInINITSel); +0317 +0318 %Here we need to figure out which original reactions (before the linear merge) +0319 %that were removed. These are all reactions with the same group ids as the removed reactions +0320 groupIdsRemoved = prepData.groupIds(ismember(prepData.refModel.rxns, deletedRxnsInINIT)); %can improve this slightly, use sel above +0321 groupIdsRemoved = groupIdsRemoved(groupIdsRemoved ~= 0);%zero means that the reaction was not grouped, all with zeros are not a group! +0322 rxnsToRem = union(prepData.refModel.rxns(ismember(prepData.groupIds,groupIdsRemoved)), deletedRxnsInINIT);%make a union here to include the ungrouped (unmerged) as well +0323 +0324 initModel = removeReactions(prepData.refModel,rxnsToRem,false,true); +0325 +0326 % remove metabolites separately to avoid removing those needed for tasks +0327 unusedMets = initModel.mets(all(initModel.S == 0,2)); +0328 initModel = removeMets(initModel, setdiff(unusedMets, prepData.essentialMetsForTasks)); +0329 +0330 %if printReport == true +0331 % printScores(initModel,'INIT model statistics',hpaData,transcrData,tissue,celltype); +0332 % printScores(removeReactions(cModel,setdiff(cModel.rxns,rxnsToRem),true,true),'Reactions deleted by INIT',hpaData,transcrData,tissue,celltype); +0333 %end +0334 +0335 %The full model has exchange reactions in it. ftINITFillGapsForAllTasks calls +0336 %ftINITFillGaps, which automatically removes exchange metabolites (because it +0337 %assumes that the reactions are constrained when appropriate). In this case the +0338 %uptakes/outputs are retrieved from the task sheet instead. To prevent +0339 %exchange reactions being used to fill gaps, they are deleted from the +0340 %reference model here. +0341 initModel.id = 'INITModel'; +0342 +0343 %If gaps in the model should be filled using a task list +0344 if ~isempty(prepData.taskStruct) +0345 %Remove exchange reactions and reactions already included in the INIT +0346 %model +0347 %We changed strategy and instead include all rxns except the exchange rxns in the ref model +0348 %But we do keep the exchange rxns that are essential. +0349 %Let's test to remove all, that should work +0350 +0351 %At this stage the model is fully connected and most of the genes with +0352 %good scores should have been included. The final gap-filling should +0353 %take the scores of the genes into account, so that "rather bad" +0354 %reactions are preferred to "very bad" reactions. However, reactions +0355 %with positive scores will be included even if they are not connected +0356 %in the current formulation. Therefore, such reactions will have to be +0357 %assigned a small negative score instead. +0358 exchRxns = getExchangeRxns(prepData.refModel); +0359 refModelNoExc = removeReactions(prepData.refModelWithBM,exchRxns,false,true); +0360 exchRxns = getExchangeRxns(initModel); +0361 initModelNoExc = removeReactions(closeModel(initModel),exchRxns,false,true); +0362 +0363 if useScoresForTasks == true +0364 %map the rxn scores to the model without exchange rxns +0365 [~,ia,ib] = intersect(refModelNoExc.rxns,prepData.refModel.rxns); +0366 rxnScores2nd = NaN(length(refModelNoExc.rxns),1); +0367 rxnScores2nd(ia) = origRxnScores(ib); +0368 %all(rxnScores2nd == refRxnScores);%should be the same, ok! +0369 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,min(rxnScores2nd,-0.1),prepData.taskStruct,paramsFT,verbose); +0370 else +0371 [outModel,addedRxnMat] = ftINITFillGapsForAllTasks(initModelNoExc,refModelNoExc,[],true,[],prepData.taskStruct,paramsFT,verbose); +0372 end +0373 %if printReport == true +0374 % printScores(outModel,'Functional model statistics',hpaData,transcrData,tissue,celltype); +0375 % printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,transcrData,tissue,celltype); +0376 %end +0377 +0378 addedRxnsForTasks = refModelNoExc.rxns(any(addedRxnMat,2)); +0379 else +0380 outModel = initModel; +0381 addedRxnMat = []; +0382 addedRxnsForTasks = {}; +0383 end +0384 +0385 % The model can now perform all the tasks defined in the task list. +0386 model = outModel; +0387 +0388 +0389 % At this stage the model will contain some exchange reactions but probably +0390 % not all (and maybe zero). This can be inconvenient, so all exchange +0391 % reactions from the reference model are added, except for those which +0392 % involve metabolites that are not in the model. +0393 +0394 %Start from the original model, and just remove the reactions that are no longer there (and keep exchange rxns). The model we got out +0395 %from the problem is not complete, it doesn't have GRPs etc. +0396 %The logic below is a bit complicated. We identify the reactions that should be removed from the full model as +0397 %reactions that have been removed in the init model except the ones that were added back. In addition, we make +0398 %sure that no exchange rxns are removed - they can be removed in the init model if they were linearly merged with other +0399 %reactions that were decided to be removed from the model. We want to keep all exchange rxns to make sure the tasks can +0400 %be performed also without manipulating the b vector in the model (which is what is done in the gap-filling). +0401 exchRxns = getExchangeRxns(prepData.refModel); +0402 deletedRxnsInINIT = setdiff(prepData.refModel.rxns,union(union(initModel.rxns, addedRxnsForTasks), exchRxns)); +0403 outModel = removeReactions(prepData.refModel, deletedRxnsInINIT, true); %we skip removing the genes for now, I'm not sure it is desirable +0404 +0405 % If requested, attempt to remove negative-score genes from the model, +0406 % depending on their role (isozyme or complex subunit) in each grRule. +0407 % See the "removeLowScoreGenes" function more more details, and to adjust +0408 % any default parameters therein. +0409 if ( removeGenes ) +0410 [~, geneScores] = scoreComplexModel(outModel,hpaData,transcrData,tissue,celltype); +0411 outModel = removeLowScoreGenes(outModel,geneScores); +0412 end +0413 +0414 +0415 model = outModel; +0416 +0417 end +0418 +0419 %This is for printing a summary of a model +0420 function [rxnS, geneS] = printScores(model,name,hpaData,transcrData,tissue,celltype) +0421 [a, b] = scoreComplexModel(model,hpaData,transcrData,tissue,celltype); +0422 rxnS = mean(a); +0423 geneS = mean(b,'omitnan'); +0424 fprintf([name ':\n']); +0425 fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']); +0426 fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']); +0427 fprintf(['\tMean gene score: ' num2str(geneS) '\n']); +0428 fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']); +0429 end +0430 +0431 function rxnsToIgnore = getRxnsFromPattern(rxnsToIgnorePattern, prepData) +0432 rxnsToIgnore = false(length(prepData.toIgnoreExch),1); +0433 if rxnsToIgnorePattern(1) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreExch; end; +0434 if rxnsToIgnorePattern(2) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreImportRxns; end; +0435 if rxnsToIgnorePattern(3) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSimpleTransp; end; +0436 if rxnsToIgnorePattern(4) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAdvTransp; end; +0437 if rxnsToIgnorePattern(5) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreSpont; end; +0438 if rxnsToIgnorePattern(6) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreS; end; +0439 if rxnsToIgnorePattern(7) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreCustomRxns; end; +0440 if rxnsToIgnorePattern(8) rxnsToIgnore = rxnsToIgnore | prepData.toIgnoreAllWithoutGPRs; end; +0441 end +0442
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/INIT/ftINITFillGapsMILP.html b/doc/INIT/ftINITFillGapsMILP.html index f2294249..888479a1 100644 --- a/doc/INIT/ftINITFillGapsMILP.html +++ b/doc/INIT/ftINITFillGapsMILP.html @@ -282,30 +282,31 @@

SOURCE CODE ^%This is weird - although it says "optimal solution found", we can get different results with different 0209 %values of the objective function, where one is more optimal than the other (pretty big difference...) -0210 %params.CSClientLog = 3;%generates a warning in gurobi, but may be of interest for other solvers -0211 -0212 % Optimize the problem -0213 res = optimizeProb(prob,params,verbose); -0214 isFeasible=checkSolution(res); -0215 -0216 if ~isFeasible -0217 x=[]; -0218 I=[]; -0219 exitFlag=-1; -0220 return; -0221 end -0222 -0223 x=res.full(1:numel(model.rxns));%the fluxes -0224 I=res.full(intsIndexes) > 10^-3;%The margin for integers in gurobi is 10^-5, not 10^-12 that was previously used! Use 10^-3 to have some margin! -0225 -0226 tmp = res.full(intsIndexes); -0227 sel = (tmp > 10^-12) & (tmp < 0.5); -0228 if sum(sel) > 0 -0229 %This may indicate that there is a problem with the tolerances in the solver -0230 disp(['ftINITFillGapsMILP: Some variables meant to be boolean in the MILP have intermediate values. Num vars: ' num2str(sum(sel))]) -0231 end -0232 -0233 end +0210 params.Threads = 1; %single-threaded Gurobi MILP for determinism, see issue #607 +0211 %params.CSClientLog = 3;%generates a warning in gurobi, but may be of interest for other solvers +0212 +0213 % Optimize the problem +0214 res = optimizeProb(prob,params,verbose); +0215 isFeasible=checkSolution(res); +0216 +0217 if ~isFeasible +0218 x=[]; +0219 I=[]; +0220 exitFlag=-1; +0221 return; +0222 end +0223 +0224 x=res.full(1:numel(model.rxns));%the fluxes +0225 I=res.full(intsIndexes) > 10^-3;%The margin for integers in gurobi is 10^-5, not 10^-12 that was previously used! Use 10^-3 to have some margin! +0226 +0227 tmp = res.full(intsIndexes); +0228 sel = (tmp > 10^-12) & (tmp < 0.5); +0229 if sum(sel) > 0 +0230 %This may indicate that there is a problem with the tolerances in the solver +0231 disp(['ftINITFillGapsMILP: Some variables meant to be boolean in the MILP have intermediate values. Num vars: ' num2str(sum(sel))]) +0232 end +0233 +0234 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/INIT/removeLowScoreGenes.html b/doc/INIT/removeLowScoreGenes.html index 34524a1e..06ff5320 100644 --- a/doc/INIT/removeLowScoreGenes.html +++ b/doc/INIT/removeLowScoreGenes.html @@ -213,153 +213,160 @@

SOURCE CODE ^% regenerate "genes" and "rxnGeneMat" model fields 0111 [genes,rxnGeneMat] = getGenesFromGrRules(newModel.grRules); -0112 newModel.genes = genes; -0113 newModel.rxnGeneMat = rxnGeneMat; -0114 -0115 % update other gene-related fields -0116 remInd = ~ismember(model.genes,newModel.genes); -0117 remGenes = model.genes(remInd); -0118 -0119 if isfield(newModel,'geneShortNames') -0120 newModel.geneShortNames(remInd) = []; -0121 end -0122 if isfield(newModel,'proteins') -0123 newModel.proteins(remInd) = []; -0124 end -0125 if isfield(newModel,'geneMiriams') -0126 newModel.geneMiriams(remInd) = []; -0127 end -0128 if isfield(newModel,'geneFrom') -0129 newModel.geneFrom(remInd) = []; -0130 end -0131 if isfield(newModel,'geneComps') -0132 newModel.geneComps(remInd) = []; -0133 end -0134 -0135 -0136 end -0137 -0138 -0139 -0140 function [updatedRule,rScore] = processSimpleRule(rule,genes,gScores,isozymeScoring,complexScoring) -0141 % Either score or modify a reaction gene rule containig only ANDs or ORs. -0142 % -0143 % If the rule contains an enzyme complex (all ANDs), the complex will be -0144 % scored based on the score of its subunits. Subunits without a score (NaN) -0145 % will be excluded from the score calculation. -0146 % -0147 % If the rule contains only isozymes (all ORs), the negative-score genes -0148 % will be removed from the rule. Isozymes without a score (NaN) will not be -0149 % removed from the rule. The resuling rule will then be scored. -0150 -0151 -0152 % get IDs and indices of genes involved in rule -0153 ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match')); -0154 [~,geneInd] = ismember(ruleGenes,genes); -0155 -0156 % rules with one or no genes remain unchanged -0157 if numel(ruleGenes) < 2 -0158 rScore = gScores(geneInd); -0159 updatedRule = rule; -0160 return -0161 end +0112 +0113 % determine which of the original genes were removed +0114 remInd = ~ismember(model.genes,genes); +0115 remGenes = model.genes(remInd); +0116 +0117 % Keep the retained genes in their original order rather than the sorted +0118 % order returned by getGenesFromGrRules. Gene removal never introduces new +0119 % genes, so the remaining genes are model.genes(~remInd). Preserving this +0120 % order ensures the gene-associated fields trimmed below (geneShortNames, +0121 % proteins, etc.), which are indexed by remInd, stay aligned with genes. +0122 newModel.genes = model.genes(~remInd); +0123 [~,reorderInd] = ismember(newModel.genes,genes); +0124 newModel.rxnGeneMat = rxnGeneMat(:,reorderInd); +0125 +0126 if isfield(newModel,'geneShortNames') +0127 newModel.geneShortNames(remInd) = []; +0128 end +0129 if isfield(newModel,'proteins') +0130 newModel.proteins(remInd) = []; +0131 end +0132 if isfield(newModel,'geneMiriams') +0133 newModel.geneMiriams(remInd) = []; +0134 end +0135 if isfield(newModel,'geneFrom') +0136 newModel.geneFrom(remInd) = []; +0137 end +0138 if isfield(newModel,'geneComps') +0139 newModel.geneComps(remInd) = []; +0140 end +0141 +0142 +0143 end +0144 +0145 +0146 +0147 function [updatedRule,rScore] = processSimpleRule(rule,genes,gScores,isozymeScoring,complexScoring) +0148 % Either score or modify a reaction gene rule containig only ANDs or ORs. +0149 % +0150 % If the rule contains an enzyme complex (all ANDs), the complex will be +0151 % scored based on the score of its subunits. Subunits without a score (NaN) +0152 % will be excluded from the score calculation. +0153 % +0154 % If the rule contains only isozymes (all ORs), the negative-score genes +0155 % will be removed from the rule. Isozymes without a score (NaN) will not be +0156 % removed from the rule. The resuling rule will then be scored. +0157 +0158 +0159 % get IDs and indices of genes involved in rule +0160 ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match')); +0161 [~,geneInd] = ismember(ruleGenes,genes); 0162 -0163 if ~contains(rule,'&') % rule contains isozymes -0164 -0165 scoreMethod = isozymeScoring; -0166 negInd = gScores(geneInd) < 0; % NaNs will return false here -0167 if all(negInd) -0168 % get the least negative gene, adding a small random value to avoid a tie -0169 [~,maxInd] = max(gScores(geneInd) + rand(size(geneInd))*(1e-8)); -0170 updatedRule = ruleGenes{maxInd}; -0171 elseif sum(~negInd) == 1 -0172 updatedRule = ruleGenes{~negInd}; -0173 else -0174 updatedRule = strjoin(ruleGenes(~negInd),' | '); -0175 if startsWith(rule,'(') -0176 updatedRule = ['(',updatedRule,')']; -0177 end -0178 end -0179 -0180 % update ruleGenes and their indices -0181 ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match')); -0182 [~,geneInd] = ismember(ruleGenes,genes); -0183 -0184 elseif ~contains(rule,'|') % rule contains enzyme complex -0185 scoreMethod = complexScoring; -0186 updatedRule = rule; -0187 else -0188 error('This function cannot handle rules with both "OR" and "AND" expressions.'); -0189 end -0190 -0191 % score rule -0192 switch lower(scoreMethod) -0193 case 'min' -0194 rScore = min(gScores(geneInd),[],'omitnan'); -0195 case 'max' -0196 rScore = max(gScores(geneInd),[],'omitnan'); -0197 case 'median' -0198 rScore = median(gScores(geneInd),'omitnan'); -0199 case 'average' -0200 rScore = mean(gScores(geneInd),'omitnan'); -0201 end -0202 -0203 end -0204 -0205 -0206 -0207 function updatedRule = processComplexRule(rule,genes,gScores,isozymeScoring,complexScoring) -0208 % Update reactions containing both AND and OR expressions. -0209 % -0210 % Negative-score genes will be removed if they are isozymic, whereas they -0211 % will not be removed if they are part of an enzyme complex. However, if -0212 % the enzyme complex has a negative score, the entire complex will be -0213 % removed, as long as it is not the only remaining element in the rule. -0214 -0215 -0216 % Specify phrases to search for in the grRule. These phrases will find -0217 % genes grouped by all ANDs (first phrase) or all ORs (second phrase). -0218 search_phrases = {'\([^&|\(\) ]+( & [^&|\(\) ]+)+\)', '\([^&|\(\) ]+( \| [^&|\(\) ]+)+\)'}; -0219 -0220 % initialize some variables -0221 subsets = {}; % subsets are groups of genes grouped by all ANDs or all ORs -0222 c = 1; % counter to keep track of the group (subset) number -0223 r_orig = rule; % record original rule to determine when it stops changing -0224 for k = 1:100 % iterate some arbitrarily high number of times -0225 for j = 1:length(search_phrases) -0226 new_subset = regexp(rule,search_phrases{j},'match')'; % extract subsets -0227 if ~isempty(new_subset) -0228 subsets = [subsets; new_subset]; % append to list of subsets -0229 subset_nums = arrayfun(@num2str,(c:length(subsets))','UniformOutput',false); % get group numbers to be assigned to the new subsets, and convert to strings -0230 rule = regexprep(rule,search_phrases{j},strcat('#',subset_nums,'#'),'once'); % replace the subsets in the expression with their group numbers (enclosed by "#"s) -0231 c = c + length(new_subset); -0232 end -0233 end -0234 if isequal(rule,r_orig) -0235 break; % stop iterating when rule stops changing -0236 else -0237 r_orig = rule; -0238 end -0239 end -0240 subsets{end+1} = rule; % add final state of rule as the last subset -0241 -0242 % score and update each subset, and append to gene list and gene scores -0243 for i = 1:numel(subsets) -0244 [subsets{i},subset_score] = processSimpleRule(subsets{i},genes,gScores,isozymeScoring,complexScoring); -0245 gScores = [gScores; subset_score]; -0246 genes = [genes; {strcat('#',num2str(i),'#')}]; -0247 end +0163 % rules with one or no genes remain unchanged +0164 if numel(ruleGenes) < 2 +0165 rScore = gScores(geneInd); +0166 updatedRule = rule; +0167 return +0168 end +0169 +0170 if ~contains(rule,'&') % rule contains isozymes +0171 +0172 scoreMethod = isozymeScoring; +0173 negInd = gScores(geneInd) < 0; % NaNs will return false here +0174 if all(negInd) +0175 % get the least negative gene, adding a small random value to avoid a tie +0176 [~,maxInd] = max(gScores(geneInd) + rand(size(geneInd))*(1e-8)); +0177 updatedRule = ruleGenes{maxInd}; +0178 elseif sum(~negInd) == 1 +0179 updatedRule = ruleGenes{~negInd}; +0180 else +0181 updatedRule = strjoin(ruleGenes(~negInd),' | '); +0182 if startsWith(rule,'(') +0183 updatedRule = ['(',updatedRule,')']; +0184 end +0185 end +0186 +0187 % update ruleGenes and their indices +0188 ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match')); +0189 [~,geneInd] = ismember(ruleGenes,genes); +0190 +0191 elseif ~contains(rule,'|') % rule contains enzyme complex +0192 scoreMethod = complexScoring; +0193 updatedRule = rule; +0194 else +0195 error('This function cannot handle rules with both "OR" and "AND" expressions.'); +0196 end +0197 +0198 % score rule +0199 switch lower(scoreMethod) +0200 case 'min' +0201 rScore = min(gScores(geneInd),[],'omitnan'); +0202 case 'max' +0203 rScore = max(gScores(geneInd),[],'omitnan'); +0204 case 'median' +0205 rScore = median(gScores(geneInd),'omitnan'); +0206 case 'average' +0207 rScore = mean(gScores(geneInd),'omitnan'); +0208 end +0209 +0210 end +0211 +0212 +0213 +0214 function updatedRule = processComplexRule(rule,genes,gScores,isozymeScoring,complexScoring) +0215 % Update reactions containing both AND and OR expressions. +0216 % +0217 % Negative-score genes will be removed if they are isozymic, whereas they +0218 % will not be removed if they are part of an enzyme complex. However, if +0219 % the enzyme complex has a negative score, the entire complex will be +0220 % removed, as long as it is not the only remaining element in the rule. +0221 +0222 +0223 % Specify phrases to search for in the grRule. These phrases will find +0224 % genes grouped by all ANDs (first phrase) or all ORs (second phrase). +0225 search_phrases = {'\([^&|\(\) ]+( & [^&|\(\) ]+)+\)', '\([^&|\(\) ]+( \| [^&|\(\) ]+)+\)'}; +0226 +0227 % initialize some variables +0228 subsets = {}; % subsets are groups of genes grouped by all ANDs or all ORs +0229 c = 1; % counter to keep track of the group (subset) number +0230 r_orig = rule; % record original rule to determine when it stops changing +0231 for k = 1:100 % iterate some arbitrarily high number of times +0232 for j = 1:length(search_phrases) +0233 new_subset = regexp(rule,search_phrases{j},'match')'; % extract subsets +0234 if ~isempty(new_subset) +0235 subsets = [subsets; new_subset]; % append to list of subsets +0236 subset_nums = arrayfun(@num2str,(c:length(subsets))','UniformOutput',false); % get group numbers to be assigned to the new subsets, and convert to strings +0237 rule = regexprep(rule,search_phrases{j},strcat('#',subset_nums,'#'),'once'); % replace the subsets in the expression with their group numbers (enclosed by "#"s) +0238 c = c + length(new_subset); +0239 end +0240 end +0241 if isequal(rule,r_orig) +0242 break; % stop iterating when rule stops changing +0243 else +0244 r_orig = rule; +0245 end +0246 end +0247 subsets{end+1} = rule; % add final state of rule as the last subset 0248 -0249 % reconstruct the rule from its updated subsets -0250 updatedRule = subsets{end}; -0251 for i = c-1:-1:1 -0252 updatedRule = regexprep(updatedRule,strcat('#',num2str(i),'#'),subsets{i}); -0253 end -0254 -0255 end -0256 -0257 -0258 +0249 % score and update each subset, and append to gene list and gene scores +0250 for i = 1:numel(subsets) +0251 [subsets{i},subset_score] = processSimpleRule(subsets{i},genes,gScores,isozymeScoring,complexScoring); +0252 gScores = [gScores; subset_score]; +0253 genes = [genes; {strcat('#',num2str(i),'#')}]; +0254 end +0255 +0256 % reconstruct the rule from its updated subsets +0257 updatedRule = subsets{end}; +0258 for i = c-1:-1:1 +0259 updatedRule = regexprep(updatedRule,strcat('#',num2str(i),'#'),subsets{i}); +0260 end +0261 +0262 end +0263 +0264 +0265
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/io/exportModel.html b/doc/io/exportModel.html index f538d90b..eacffe43 100644 --- a/doc/io/exportModel.html +++ b/doc/io/exportModel.html @@ -627,239 +627,242 @@

SOURCE CODE ^'UniformOutput', false); 0565 model.subSystems(cellfun(@isempty, model.subSystems)) = {{}}; 0566 -0567 % If some entries contain string scalars by mistake, coerce them: -0568 model.subSystems = cellfun(@(c) cellfun(@char, c, 'UniformOutput', false), model.subSystems, 'UniformOutput', false); -0569 -0570 % === 2) Flatten once: names and their reaction indices (vectorized) === -0571 flatNames = vertcat(model.subSystems{:}); % 1×M cellstr of all subsystem labels -0572 if isempty(flatNames) -0573 % Nothing to do: no subsystems present -0574 return -0575 end -0576 -0577 counts = cellfun(@numel, model.subSystems); % reactions -> how many subsystems -0578 % For each reaction r, repeat r exactly counts(r) times -0579 flatIdx = arrayfun(@(r,c) repmat(r, c, 1), (1:numel(model.subSystems)).', counts, 'UniformOutput', false); -0580 flatIdx = vertcat(flatIdx{:}); % M×1 vector of reaction indices -0581 -0582 % === 3) Group in one shot: unique subsystems + members per subsystem === -0583 [subSystems, ~, g] = unique(flatNames, 'stable'); % stable preserves first appearance -0584 membersIdx = accumarray(g, flatIdx, [], @(v){v}); % cell: one vector of r-idx per group -0585 nSubs = numel(subSystems); -0586 -0587 % === 4) Preallocate and build SBML groups (single simple loop) === -0588 modelSBML.groups_group(1:nSubs) = grpTemplate; -0589 groupIDs = "group" + (1:nSubs); -0590 -0591 if isfield(grpTemplate,'groups_member') -0592 memTemplate = grpTemplate.groups_member; -0593 else -0594 memTemplate = struct('groups_idRef',''); -0595 end -0596 -0597 for i = 1:nSubs -0598 rIdx = membersIdx{i}; -0599 groupRXNs = model.rxns(rIdx); -0600 cgroup = grpTemplate; -0601 -0602 % Preallocate and fill members -0603 nM = numel(groupRXNs); -0604 if nM > 0 -0605 cgroup.groups_member(1:nM) = memTemplate; -0606 [cgroup.groups_member.groups_idRef] = deal(groupRXNs{:}); -0607 else -0608 cgroup.groups_member = repmat(memTemplate, 0, 1); -0609 end -0610 -0611 cgroup.groups_id = char(groupIDs(i)); % keep as char for SBML compatibility -0612 cgroup.groups_name = subSystems{i}; -0613 modelSBML.groups_group(i) = cgroup; -0614 end -0615 end -0616 -0617 -0618 %Prepare fbc_objective subfield +0567 % If some entries contain string scalars by mistake, coerce them. Also +0568 % force each entry to a column cell array, so that reactions with +0569 % different numbers of subsystems can be concatenated below (otherwise +0570 % mixing e.g. 1x1 and 1x2 entries breaks vertcat): +0571 model.subSystems = cellfun(@(c) reshape(cellfun(@char, c, 'UniformOutput', false), [], 1), model.subSystems, 'UniformOutput', false); +0572 +0573 % === 2) Flatten once: names and their reaction indices (vectorized) === +0574 flatNames = vertcat(model.subSystems{:}); % Mx1 cellstr of all subsystem labels +0575 if isempty(flatNames) +0576 % Nothing to do: no subsystems present +0577 return +0578 end +0579 +0580 counts = cellfun(@numel, model.subSystems); % reactions -> how many subsystems +0581 % For each reaction r, repeat r exactly counts(r) times +0582 flatIdx = arrayfun(@(r,c) repmat(r, c, 1), (1:numel(model.subSystems)).', counts, 'UniformOutput', false); +0583 flatIdx = vertcat(flatIdx{:}); % M×1 vector of reaction indices +0584 +0585 % === 3) Group in one shot: unique subsystems + members per subsystem === +0586 [subSystems, ~, g] = unique(flatNames, 'stable'); % stable preserves first appearance +0587 membersIdx = accumarray(g, flatIdx, [], @(v){v}); % cell: one vector of r-idx per group +0588 nSubs = numel(subSystems); +0589 +0590 % === 4) Preallocate and build SBML groups (single simple loop) === +0591 modelSBML.groups_group(1:nSubs) = grpTemplate; +0592 groupIDs = "group" + (1:nSubs); +0593 +0594 if isfield(grpTemplate,'groups_member') +0595 memTemplate = grpTemplate.groups_member; +0596 else +0597 memTemplate = struct('groups_idRef',''); +0598 end +0599 +0600 for i = 1:nSubs +0601 rIdx = membersIdx{i}; +0602 groupRXNs = model.rxns(rIdx); +0603 cgroup = grpTemplate; +0604 +0605 % Preallocate and fill members +0606 nM = numel(groupRXNs); +0607 if nM > 0 +0608 cgroup.groups_member(1:nM) = memTemplate; +0609 [cgroup.groups_member.groups_idRef] = deal(groupRXNs{:}); +0610 else +0611 cgroup.groups_member = repmat(memTemplate, 0, 1); +0612 end +0613 +0614 cgroup.groups_id = char(groupIDs(i)); % keep as char for SBML compatibility +0615 cgroup.groups_name = subSystems{i}; +0616 modelSBML.groups_group(i) = cgroup; +0617 end +0618 end 0619 -0620 modelSBML.fbc_objective.fbc_type='maximize'; -0621 modelSBML.fbc_objective.fbc_id='obj'; +0620 +0621 %Prepare fbc_objective subfield 0622 -0623 ind=find(model.c); -0624 -0625 if isempty(ind) -0626 modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0; -0627 else -0628 for i=1:length(ind) -0629 %Copy the default values to the next index as long as it is not the -0630 %last one -0631 if i<numel(ind) -0632 modelSBML.reaction(i+1)=modelSBML.reaction(i); -0633 end -0634 values=model.c(model.c~=0); -0635 modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id; -0636 modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i); -0637 modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1; -0638 end -0639 end -0640 -0641 modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id; -0642 -0643 fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))]; -0644 if modelHasSubsystems -0645 groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))]; -0646 modelSBML.namespaces=struct('prefix',{'','fbc','groups'},... -0647 'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],... -0648 fbcStr,groupStr}); -0649 else -0650 modelSBML.namespaces=struct('prefix',{'','fbc'},... -0651 'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],... -0652 fbcStr}); -0653 end -0654 -0655 if sbmlPackageVersions(1) == 2 -0656 modelSBML.fbc_strict=1; -0657 modelSBML.isSetfbc_strict = 1; -0658 end -0659 -0660 modelSBML.rule=[]; -0661 modelSBML.constraint=[]; +0623 modelSBML.fbc_objective.fbc_type='maximize'; +0624 modelSBML.fbc_objective.fbc_id='obj'; +0625 +0626 ind=find(model.c); +0627 +0628 if isempty(ind) +0629 modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0; +0630 else +0631 for i=1:length(ind) +0632 %Copy the default values to the next index as long as it is not the +0633 %last one +0634 if i<numel(ind) +0635 modelSBML.reaction(i+1)=modelSBML.reaction(i); +0636 end +0637 values=model.c(model.c~=0); +0638 modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id; +0639 modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i); +0640 modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1; +0641 end +0642 end +0643 +0644 modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id; +0645 +0646 fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))]; +0647 if modelHasSubsystems +0648 groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))]; +0649 modelSBML.namespaces=struct('prefix',{'','fbc','groups'},... +0650 'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],... +0651 fbcStr,groupStr}); +0652 else +0653 modelSBML.namespaces=struct('prefix',{'','fbc'},... +0654 'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],... +0655 fbcStr}); +0656 end +0657 +0658 if sbmlPackageVersions(1) == 2 +0659 modelSBML.fbc_strict=1; +0660 modelSBML.isSetfbc_strict = 1; +0661 end 0662 -0663 [ravenDir,prevDir]=findRAVENroot(); -0664 fileName=checkFileExistence(fileName,1,true,false); +0663 modelSBML.rule=[]; +0664 modelSBML.constraint=[]; 0665 -0666 OutputSBML_RAVEN(modelSBML,fileName,1,0,[1,0]); -0667 end +0666 [ravenDir,prevDir]=findRAVENroot(); +0667 fileName=checkFileExistence(fileName,1,true,false); 0668 -0669 -0670 function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions) -0671 %Returns the blank SBML model structure by using appropriate libSBML -0672 %functions. This creates structure by considering three levels -0673 -0674 sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0675 sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0669 OutputSBML_RAVEN(modelSBML,fileName,1,0,[1,0]); +0670 end +0671 +0672 +0673 function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions) +0674 %Returns the blank SBML model structure by using appropriate libSBML +0675 %functions. This creates structure by considering three levels 0676 -0677 for i=1:numel(sbmlFieldNames) -0678 modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i}; -0679 sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0680 sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0681 if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment') -0682 for j=1:numel(sbmlSubfieldNames) -0683 modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j}; -0684 sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0685 sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0686 if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw') -0687 for k=1:numel(sbmlSubsubfieldNames) -0688 %'compartment' and 'species' fields are not supposed to -0689 %have their standalone structures if they are subfields -0690 %or subsubfields -0691 if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species') -0692 modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k}; -0693 end -0694 %If it is fbc_association in the third level, we need -0695 %to establish the fourth level, since libSBML requires -0696 %it -0697 if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association') -0698 fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0699 fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0700 for l=1:numel(fbc_associationFieldNames) -0701 modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l}; -0702 end -0703 end -0704 end -0705 end -0706 end -0707 end -0708 if ~isstruct(modelSBML.(sbmlFieldNames{1,i})) -0709 modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i}; +0677 sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0678 sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0679 +0680 for i=1:numel(sbmlFieldNames) +0681 modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i}; +0682 sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0683 sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0684 if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment') +0685 for j=1:numel(sbmlSubfieldNames) +0686 modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j}; +0687 sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0688 sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0689 if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw') +0690 for k=1:numel(sbmlSubsubfieldNames) +0691 %'compartment' and 'species' fields are not supposed to +0692 %have their standalone structures if they are subfields +0693 %or subsubfields +0694 if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species') +0695 modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k}; +0696 end +0697 %If it is fbc_association in the third level, we need +0698 %to establish the fourth level, since libSBML requires +0699 %it +0700 if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association') +0701 fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0702 fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0703 for l=1:numel(fbc_associationFieldNames) +0704 modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l}; +0705 end +0706 end +0707 end +0708 end +0709 end 0710 end -0711 end -0712 -0713 modelSBML.unitDefinition.id='mmol_per_gDW_per_hr'; -0714 -0715 unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); -0716 unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0711 if ~isstruct(modelSBML.(sbmlFieldNames{1,i})) +0712 modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i}; +0713 end +0714 end +0715 +0716 modelSBML.unitDefinition.id='mmol_per_gDW_per_hr'; 0717 -0718 kinds={'mole','gram','second'}; -0719 exponents=[1 -1 -1]; -0720 scales=[-3 0 0]; -0721 multipliers=[1 1 1*60*60]; -0722 -0723 for i=1:numel(unitFieldNames) -0724 modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i}; -0725 for j=1:3 -0726 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i}; -0727 if strcmp(unitFieldNames{1,i},'kind') -0728 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j}; -0729 elseif strcmp(unitFieldNames{1,i},'exponent') -0730 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j); -0731 elseif strcmp(unitFieldNames{1,i},'scale') -0732 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j); -0733 elseif strcmp(unitFieldNames{1,i},'multiplier') -0734 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j); -0735 end -0736 end -0737 end -0738 end -0739 -0740 function miriamString=getMiriam(miriamStruct) -0741 %Returns a string with list elements for a miriam structure ('<rdf:li -0742 %rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This -0743 %is just to speed up things since this is done many times during the -0744 %exporting -0745 -0746 miriamString=''; -0747 if isfield(miriamStruct,'name') -0748 for i=1:numel(miriamStruct.name) -0749 miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>']; -0750 end -0751 end -0752 end -0753 -0754 function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i) -0755 %This function provides reactants and products for particular reaction. The -0756 %function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679 -0757 -0758 met_idx = find(model.S(:, i)); -0759 tmp_Rxn.product=[]; -0760 tmp_Rxn.reactant=[]; -0761 for j_met=1:size(met_idx,1) -0762 tmp_idx = met_idx(j_met,1); -0763 sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id; -0764 met_stoich = model.S(tmp_idx, i); -0765 sbml_tmp_species_ref.stoichiometry = abs(met_stoich); -0766 sbml_tmp_species_ref.isSetStoichiometry=1; -0767 sbml_tmp_species_ref.constant=1; -0768 if (met_stoich > 0) -0769 tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ]; -0770 else -0771 tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref]; -0772 end -0773 end -0774 end -0775 -0776 function vecT = columnVector(vec) -0777 % Code below taken from COBRA Toolbox under GNU General Public License v3.0 -0778 % license file in readme/GPL.MD. -0779 % -0780 % Converts a vector to a column vector -0781 % -0782 % USAGE: -0783 % -0784 % vecT = columnVector(vec) -0785 % -0786 % INPUT: -0787 % vec: a vector +0718 unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0719 unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions); +0720 +0721 kinds={'mole','gram','second'}; +0722 exponents=[1 -1 -1]; +0723 scales=[-3 0 0]; +0724 multipliers=[1 1 1*60*60]; +0725 +0726 for i=1:numel(unitFieldNames) +0727 modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i}; +0728 for j=1:3 +0729 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i}; +0730 if strcmp(unitFieldNames{1,i},'kind') +0731 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j}; +0732 elseif strcmp(unitFieldNames{1,i},'exponent') +0733 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j); +0734 elseif strcmp(unitFieldNames{1,i},'scale') +0735 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j); +0736 elseif strcmp(unitFieldNames{1,i},'multiplier') +0737 modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j); +0738 end +0739 end +0740 end +0741 end +0742 +0743 function miriamString=getMiriam(miriamStruct) +0744 %Returns a string with list elements for a miriam structure ('<rdf:li +0745 %rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This +0746 %is just to speed up things since this is done many times during the +0747 %exporting +0748 +0749 miriamString=''; +0750 if isfield(miriamStruct,'name') +0751 for i=1:numel(miriamStruct.name) +0752 miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>']; +0753 end +0754 end +0755 end +0756 +0757 function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i) +0758 %This function provides reactants and products for particular reaction. The +0759 %function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679 +0760 +0761 met_idx = find(model.S(:, i)); +0762 tmp_Rxn.product=[]; +0763 tmp_Rxn.reactant=[]; +0764 for j_met=1:size(met_idx,1) +0765 tmp_idx = met_idx(j_met,1); +0766 sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id; +0767 met_stoich = model.S(tmp_idx, i); +0768 sbml_tmp_species_ref.stoichiometry = abs(met_stoich); +0769 sbml_tmp_species_ref.isSetStoichiometry=1; +0770 sbml_tmp_species_ref.constant=1; +0771 if (met_stoich > 0) +0772 tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ]; +0773 else +0774 tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref]; +0775 end +0776 end +0777 end +0778 +0779 function vecT = columnVector(vec) +0780 % Code below taken from COBRA Toolbox under GNU General Public License v3.0 +0781 % license file in readme/GPL.MD. +0782 % +0783 % Converts a vector to a column vector +0784 % +0785 % USAGE: +0786 % +0787 % vecT = columnVector(vec) 0788 % -0789 % OUTPUT: -0790 % vecT: a column vector -0791 -0792 [n, m] = size(vec); -0793 -0794 if n < m -0795 vecT = vec'; -0796 else -0797 vecT = vec; -0798 end -0799 end +0789 % INPUT: +0790 % vec: a vector +0791 % +0792 % OUTPUT: +0793 % vecT: a column vector +0794 +0795 [n, m] = size(vec); +0796 +0797 if n < m +0798 vecT = vec'; +0799 else +0800 vecT = vec; +0801 end +0802 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/testing/unit_tests/importExportTests.html b/doc/testing/unit_tests/importExportTests.html index 60dc50cf..00201adb 100644 --- a/doc/testing/unit_tests/importExportTests.html +++ b/doc/testing/unit_tests/importExportTests.html @@ -42,7 +42,7 @@

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

+
  • function testExcelImport(testCase)
  • function testSBMLImport(testCase)
  • function testYAMLimport(testCase)
  • function testExcelExport(testCase)
  • function testSBMLExport(testCase)
  • function testSBMLExportNestedSubSystems(testCase)
  • function testYAMLexport(testCase)

    • SOURCE CODE ^

    0001 %run this test case with the command
@@ -106,18 +106,48 @@ 
    SOURCE CODE ^'testing','unit_tests','test_data','_test.xml'));
 0060 end
 0061 
-0062 function testYAMLexport(testCase)
-0063 sourceDir=fileparts(fileparts(fileparts(which(mfilename))));
-0064 load(fullfile(sourceDir,'tutorial','empty.mat'), 'emptyModel');
-0065 evalc('writeYAMLmodel(emptyModel,fullfile(sourceDir,''testing'',''unit_tests'',''test_data'',''_test.yml''))');
-0066 %File will not be exactly equal as it contains the current date and time,
-0067 %so md5 or similar would not work. Just check whether file is reasonably
-0068 %sized.
-0069 s = dir(fullfile(sourceDir,'testing','unit_tests','test_data','_test.yml'));         
-0070 filesize = s.bytes;
-0071 verifyTrue(testCase,filesize>1290);
-0072 delete(fullfile(sourceDir,'testing','unit_tests','test_data','_test.yml'));
-0073 end
    +0062 function testSBMLExportNestedSubSystems(testCase) +0063 %Regression test: a model where reactions have differing numbers of +0064 %subsystems (nested cell-of-cells, e.g. one reaction in two subsystems and +0065 %others in one) must be exportable to SBML. Previously exportModel failed +0066 %with "Dimensions of arrays being concatenated are not consistent" because +0067 %the subSystems flattening could not concatenate entries of unequal length. +0068 sourceDir=fileparts(fileparts(fileparts(which(mfilename)))); +0069 load(fullfile(sourceDir,'testing','unit_tests','test_data','ecoli_textbook.mat'), 'model'); +0070 +0071 %Assign nested subSystems: most reactions get a single subsystem, but a few +0072 %get multiple, which is the scenario that triggered the bug. +0073 nRxns=numel(model.rxns); +0074 model.subSystems=repmat({{'Subsystem A'}},nRxns,1); +0075 model.subSystems{1}={'Subsystem A','Subsystem B'}; +0076 model.subSystems{2}={'Subsystem B','Subsystem C'}; +0077 +0078 tmpFile=fullfile(sourceDir,'testing','unit_tests','test_data','_testNested.xml'); +0079 %This call previously errored; verify it completes and round-trips. +0080 evalc('exportModel(model,tmpFile)'); +0081 evalc('modelImported=importModel(tmpFile)'); +0082 delete(tmpFile); +0083 +0084 %Subsystems should be preserved (order within a reaction is not guaranteed) +0085 srt=@(c) sort(reshape(c,1,[])); +0086 for i=1:3 +0087 j=find(strcmp(modelImported.rxns,model.rxns{i})); +0088 verifyEqual(testCase,srt(modelImported.subSystems{j}),srt(model.subSystems{i})); +0089 end +0090 end +0091 +0092 function testYAMLexport(testCase) +0093 sourceDir=fileparts(fileparts(fileparts(which(mfilename)))); +0094 load(fullfile(sourceDir,'tutorial','empty.mat'), 'emptyModel'); +0095 evalc('writeYAMLmodel(emptyModel,fullfile(sourceDir,''testing'',''unit_tests'',''test_data'',''_test.yml''))'); +0096 %File will not be exactly equal as it contains the current date and time, +0097 %so md5 or similar would not work. Just check whether file is reasonably +0098 %sized. +0099 s = dir(fullfile(sourceDir,'testing','unit_tests','test_data','_test.yml')); +0100 filesize = s.bytes; +0101 verifyTrue(testCase,filesize>1290); +0102 delete(fullfile(sourceDir,'testing','unit_tests','test_data','_test.yml')); +0103 end
    Generated by m2html © 2005
    \ No newline at end of file