diff --git a/INIT/INITStepDesc.m b/INIT/INITStepDesc.m index 91ec541c..22e8a67d 100644 --- a/INIT/INITStepDesc.m +++ b/INIT/INITStepDesc.m @@ -24,7 +24,9 @@ % .mets Names of metabolites to remove % .compsToKeep Compartments for which metabolites should be kept. MILPParams %Cell array of MILPparams - dictates how many iterations that will be run in this step. - %Typically, MIPGap and TimeLimit is specified + %Typically, MIPGap and TimeLimit is specified. ftINIT additionally + %defaults Threads to 1 (single-threaded Gurobi) for determinism; + %set Threads here to override (0 = use all cores). AbsMIPGaps %If the objective is close to zero, a percentage of that is very small. %Therefore, also set an absolut value for this (typically 10 or 20). %For practical reasons, the first number is not used diff --git a/INIT/ftINIT.m b/INIT/ftINIT.m index 29654b29..ab5f3e03 100644 --- a/INIT/ftINIT.m +++ b/INIT/ftINIT.m @@ -233,6 +233,13 @@ if ~isfield(params, 'TimeLimit') params.TimeLimit = 5000; end + + %Default to single-threaded Gurobi MILP solving. Multi-threaded Gurobi + %can non-deterministically report the MILP as infeasible (issue #607). + %Override by setting 'Threads' in the step's MILPParams (0 = all cores). + if ~isfield(params, 'Threads') + params.Threads = 1; + end if ~first %There is sometimes a problem with that the objective function becomes close to zero, diff --git a/INIT/ftINITFillGapsMILP.m b/INIT/ftINITFillGapsMILP.m index cc66f33c..3eb51ac7 100644 --- a/INIT/ftINITFillGapsMILP.m +++ b/INIT/ftINITFillGapsMILP.m @@ -207,6 +207,7 @@ params.TimeLimit = 300; params.Seed = 26;%This is weird - although it says "optimal solution found", we can get different results with different %values of the objective function, where one is more optimal than the other (pretty big difference...) +params.Threads = 1; %single-threaded Gurobi MILP for determinism, see issue #607 %params.CSClientLog = 3;%generates a warning in gurobi, but may be of interest for other solvers % Optimize the problem 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 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 diff --git a/io/exportModel.m b/io/exportModel.m index 27c9bb6f..5deeea08 100755 --- a/io/exportModel.m +++ b/io/exportModel.m @@ -564,11 +564,14 @@ function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds) model.subSystems(isChar) = cellfun(@(s){s}, model.subSystems(isChar), 'UniformOutput', false); model.subSystems(cellfun(@isempty, model.subSystems)) = {{}}; - % If some entries contain string scalars by mistake, coerce them: - model.subSystems = cellfun(@(c) cellfun(@char, c, 'UniformOutput', false), model.subSystems, 'UniformOutput', false); + % If some entries contain string scalars by mistake, coerce them. Also + % force each entry to a column cell array, so that reactions with + % different numbers of subsystems can be concatenated below (otherwise + % mixing e.g. 1x1 and 1x2 entries breaks vertcat): + model.subSystems = cellfun(@(c) reshape(cellfun(@char, c, 'UniformOutput', false), [], 1), model.subSystems, 'UniformOutput', false); % === 2) Flatten once: names and their reaction indices (vectorized) === - flatNames = vertcat(model.subSystems{:}); % 1×M cellstr of all subsystem labels + flatNames = vertcat(model.subSystems{:}); % Mx1 cellstr of all subsystem labels if isempty(flatNames) % Nothing to do: no subsystems present return diff --git a/testing/unit_tests/importExportTests.m b/testing/unit_tests/importExportTests.m index ff670f57..aadaf82b 100755 --- a/testing/unit_tests/importExportTests.m +++ b/testing/unit_tests/importExportTests.m @@ -59,6 +59,36 @@ function testSBMLExport(testCase) delete(fullfile(sourceDir,'testing','unit_tests','test_data','_test.xml')); end +function testSBMLExportNestedSubSystems(testCase) +%Regression test: a model where reactions have differing numbers of +%subsystems (nested cell-of-cells, e.g. one reaction in two subsystems and +%others in one) must be exportable to SBML. Previously exportModel failed +%with "Dimensions of arrays being concatenated are not consistent" because +%the subSystems flattening could not concatenate entries of unequal length. +sourceDir=fileparts(fileparts(fileparts(which(mfilename)))); +load(fullfile(sourceDir,'testing','unit_tests','test_data','ecoli_textbook.mat'), 'model'); + +%Assign nested subSystems: most reactions get a single subsystem, but a few +%get multiple, which is the scenario that triggered the bug. +nRxns=numel(model.rxns); +model.subSystems=repmat({{'Subsystem A'}},nRxns,1); +model.subSystems{1}={'Subsystem A','Subsystem B'}; +model.subSystems{2}={'Subsystem B','Subsystem C'}; + +tmpFile=fullfile(sourceDir,'testing','unit_tests','test_data','_testNested.xml'); +%This call previously errored; verify it completes and round-trips. +evalc('exportModel(model,tmpFile)'); +evalc('modelImported=importModel(tmpFile)'); +delete(tmpFile); + +%Subsystems should be preserved (order within a reaction is not guaranteed) +srt=@(c) sort(reshape(c,1,[])); +for i=1:3 + j=find(strcmp(modelImported.rxns,model.rxns{i})); + verifyEqual(testCase,srt(modelImported.subSystems{j}),srt(model.subSystems{i})); +end +end + function testYAMLexport(testCase) sourceDir=fileparts(fileparts(fileparts(which(mfilename)))); load(fullfile(sourceDir,'tutorial','empty.mat'), 'emptyModel');