From da342345318ae1a577e9461f5be6b7ac53072462 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 17:01:12 +0100 Subject: [PATCH] fix: four bugs found while rewriting the help blocks, corrected examples and comments - `renameModelGenes` reads the gene table from a file path (the file name was passed to `%s` as a table). - `getIndexes` with type `metnames` and `returnLogical` true raises `RAVEN:badInput` (each name can match several metabolites). - `checkProduction` returns a model with excretion reactions only when `checkNeededForProduction` is true. - `findGeneDeletions` returns `grRatioMuts` with one entry per gene pair for double deletions. - `checkRaven` names the offline bundle `RAVEN--full.zip`. The examples of `gapReport`, `fitTasks`, `printFluxes` and `predictLocalization` use the documented arguments. The comments of `simplifyModel`, `setExchangeBounds`, `closeModel` and `ftINIT` state the current behaviour, and the help of `getModelFromKEGG` states which reactions the keep options retain. Each fix has a function test. --- CHANGELOG.md | 25 ++++++++++++++++++++++ checkRaven.m | 2 +- src/INIT/ftINIT.m | 8 +++---- src/analysis/findGeneDeletions.m | 2 +- src/curation/renameModelGenes.m | 9 ++++---- src/gapfilling/checkProduction.m | 13 ++++++----- src/gapfilling/fitTasks.m | 4 ++-- src/gapfilling/gapReport.m | 5 +++-- src/localization/predictLocalization.m | 3 ++- src/manipulation/closeModel.m | 14 ++++-------- src/manipulation/setExchangeBounds.m | 6 ++---- src/manipulation/simplifyModel.m | 7 +++--- src/queries/getIndexes.m | 9 +++++++- src/queries/printFluxes.m | 4 ++-- src/reconstruction/kegg/getModelFromKEGG.m | 4 +++- tests/function_tests/tAnalysis.m | 18 ++++++++++++++++ tests/function_tests/tGapfilling.m | 22 +++++++++++++++++++ tests/function_tests/tGenomeData.m | 19 ++++++++++++++++ tests/function_tests/tQueries.m | 13 +++++++++++ 19 files changed, 146 insertions(+), 41 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cb5d234a..3f3dddaf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,27 @@ * `gapFillTopological`: without `targets` and without an objective reaction it raises `RAVEN:badInput` (was: an error without identifier). Help states the default seeds and targets. +* `checkProduction` with `checkNeededForProduction`: the returned `model` holds excretion + reactions only (was: the excretion reaction of the last metabolite in `notProduced` was + returned as an uptake reaction). This matters for code that uses the fifth output. + +### Curation + +* `renameModelGenes` with a file path as `geneTable`: the table is read and the genes are + renamed (was: an error `MATLAB:invalidConversion` when the message that names the file + was printed). + +### Queries + +* `getIndexes` with type `'metnames'` and `returnLogical` true: raises `RAVEN:badInput` (was: + the error "Unable to use a value of type cell as an index"). A name can match several + metabolites, so no single logical vector is defined. + +### Analysis + +* `findGeneDeletions` with `'dgd'`: `grRatioMuts` has one element per gene pair, as the help + states (was: twice as many, the second half zero). This matters for code that indexes + `grRatioMuts` by the number of its elements. ### Context-specific models (ftINIT / tINIT) @@ -57,6 +78,10 @@ case, `&`, `&&`, `|`, `||`, operators as whole words, so `G1&G2` is one gene id), and that a malformed rule raises `RAVEN:badGrRule`. +* `checkRaven`: when the binaries cannot be downloaded, the message names the + `RAVEN--full.zip` bundle attached to each release for an offline install (was: + "the `*-binaries` RAVEN release"). + ### Documentation * Help blocks of 176 files are rewritten in impersonal, present-tense wording with British spelling; the code is unchanged. Statements that differed from the code are corrected, for example `simplifyModel` `deleteDuplicates` (first occurrence kept), `removeBadRxns` `refModel`, `replaceMets` outputs, `getMinNrFluxes` `params`, `checkRaven` `developMode` and the `downloadRavenBinaries` offline bundle. `reverseRxns` has the RAVEN 3 help structure. diff --git a/checkRaven.m b/checkRaven.m index 380c1e2f..d207bec2 100755 --- a/checkRaven.m +++ b/checkRaven.m @@ -257,7 +257,7 @@ downloadRavenBinaries({'blast+','diamond','hmmer'}); catch ME printOrange([' Could not download binaries from toolbox-assets: ' ME.message '\n' ... - ' (offline? fetch the "*-binaries" RAVEN release for an offline install)\n']); + ' (offline? install from the RAVEN--full.zip bundle attached to each RAVEN release)\n']); end % Run each bundled binary with a version/help flag to confirm it % executes on this system, rather than running the full test cases. diff --git a/src/INIT/ftINIT.m b/src/INIT/ftINIT.m index 6b5ba6e4..a81496c4 100644 --- a/src/INIT/ftINIT.m +++ b/src/INIT/ftINIT.m @@ -285,10 +285,10 @@ lastError = []; if ~isempty(proveAbsGap) %One solve proven to this fixed absolute gap, replacing the relative-gap - %escalation below. The escalation's final near-zero-objective run otherwise - %accepts an arbitrary within-gap incumbent -- measured on genome-scale - %Human-GEM to cost 2.0-4.0 objective units (raven-toolbox's prove_abs_gap, - %raven-gecko-parity#104). Pairs naturally with resolveTies. + %escalation below. The final near-zero-objective run of the escalation + %accepts an arbitrary incumbent within the gap; on genome-scale Human-GEM + %the step 1 objective of the escalation ended 4.0 to 12.0 units below the + %proven optimum (raven-toolbox's prove_abs_gap, raven-gecko-parity#104). if ~isempty(stp.MILPParams) && isfield(stp.MILPParams{1}, 'TimeLimit') params = struct('TimeLimit', stp.MILPParams{1}.TimeLimit); else diff --git a/src/analysis/findGeneDeletions.m b/src/analysis/findGeneDeletions.m index 9ad6d278..e944e2de 100755 --- a/src/analysis/findGeneDeletions.m +++ b/src/analysis/findGeneDeletions.m @@ -113,7 +113,7 @@ [~, I]=ismember(originalGenes(details==1),model.genes); genesToModify=nchoosek(I,2); genes=geneMapping(genesToModify); - grRatioMuts=zeros(1,numel(genes)); + grRatioMuts=zeros(1,size(genesToModify,1)); fluxes=sparse(numel(model.rxns),size(genesToModify,1)); for i=1:size(genesToModify,1) tempModel=removeGenes(model,genesToModify(i,:),false,false,false); diff --git a/src/curation/renameModelGenes.m b/src/curation/renameModelGenes.m index 9dc2a32f..12425541 100644 --- a/src/curation/renameModelGenes.m +++ b/src/curation/renameModelGenes.m @@ -80,14 +80,15 @@ % Load table if a file path was supplied if ischar(geneTable) || isstring(geneTable) - if ~isfile(geneTable) + tableFile = char(geneTable); + if ~isfile(tableFile) error('renameModelGenes:fileNotFound', ... - 'Table file not found: %s', geneTable); + 'Table file not found: %s', tableFile); end - geneTable = readtable(geneTable, 'Delimiter', '\t', 'FileType', 'text', ... + geneTable = readtable(tableFile, 'Delimiter', '\t', 'FileType', 'text', ... 'TextType', 'char'); fprintf('Loaded gene table from file: %s (%d rows)\n', ... - geneTable, height(geneTable)); + tableFile, height(geneTable)); end % Validate column names diff --git a/src/gapfilling/checkProduction.m b/src/gapfilling/checkProduction.m index 157bede6..0f5cec0c 100755 --- a/src/gapfilling/checkProduction.m +++ b/src/gapfilling/checkProduction.m @@ -52,11 +52,8 @@ % checkNeededForProduction is false. % model : struct % model structure with an excretion reaction added for each metabolite -% in the compartments of excretionFromCompartments. When -% checkNeededForProduction is true, the excretion reaction of the last -% metabolite in notProduced is returned as an uptake reaction, because -% the loop that flips each reaction into an uptake reaction does not -% reset the last one. +% in the compartments of excretionFromCompartments. All added reactions +% are excretion reactions, also when checkNeededForProduction is true. % % Examples % -------- @@ -104,6 +101,12 @@ %be produced now neededForProductionMat(i,:)=haveFlux(model,10^-5,addedRxns(notProduced)); end + if ~isempty(notProduced) + %Reset the last iteration, so that the returned model has excretion + %reactions only + lastCol=numel(model.rxns)-numel(addedRxns)+notProduced(end); + model.S(:,lastCol)=model.S(:,lastCol)*-1; + end %Calculate the smallest number of metabolites that must be connected to %make everything connected and return their names diff --git a/src/gapfilling/fitTasks.m b/src/gapfilling/fitTasks.m index b45beb94..6ea251a3 100755 --- a/src/gapfilling/fitTasks.m +++ b/src/gapfilling/fitTasks.m @@ -65,8 +65,8 @@ % % Examples % -------- -% [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,... -% rxnScores,taskStructure,params); +% [outModel, addedRxns, failedTasks]=fitTasks(model,refModel,inputFile,... +% 'printOutput',false,'rxnScores',rxnScores,'params',params); p=parseRAVENargs(varargin, {'printOutput',true; ... 'rxnScores',[]; ... diff --git a/src/gapfilling/gapReport.m b/src/gapfilling/gapReport.m index f8a006d4..c5b1d1c8 100755 --- a/src/gapfilling/gapReport.m +++ b/src/gapfilling/gapReport.m @@ -69,8 +69,9 @@ % Examples % -------- % [noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, ... -% minToConnect, neededForProductionMat, connectedFromTemplates, ... -% addedFromTemplates] = gapReport(model, templateModels); +% minToConnect, neededForProductionMat, canProduceWithoutInput, ... +% canConsumeWithoutOutput, connectedFromTemplates, addedFromTemplates] = ... +% gapReport(model, 'templateModels', templateModels); p=parseRAVENargs(varargin, {'templateModels',[]}); templateModels=p.templateModels; diff --git a/src/localization/predictLocalization.m b/src/localization/predictLocalization.m index 05c1ca12..015360e2 100755 --- a/src/localization/predictLocalization.m +++ b/src/localization/predictLocalization.m @@ -64,7 +64,8 @@ % -------- % [outModel, geneLocalization, transportStruct, scores, removedRxns] = ... % predictLocalization(model, GSS, defaultCompartment, ... -% transportCost, maxTime, plotResults); +% 'transportCost', transportCost, 'maxTime', maxTime, ... +% 'plotResults', plotResults); p=parseRAVENargs(varargin, {'transportCost',[]; 'maxTime',[]; 'plotResults',[]}); transportCost=p.transportCost; diff --git a/src/manipulation/closeModel.m b/src/manipulation/closeModel.m index 5def2209..d3bea092 100755 --- a/src/manipulation/closeModel.m +++ b/src/manipulation/closeModel.m @@ -36,15 +36,9 @@ col=closedModel.S(:,i); %A boundary/exchange/sink/demand reaction has metabolites on only one %side, same definition as getExchangeRxns: no substrates or no - %products. The previous coefficient-magnitude rule (all nonzero - %coefficients summing to 1 in absolute value) missed a scaled - %single-metabolite sink such as "2 A =>" (sum 2) outright, and only - %matched a multi-metabolite one-sided reaction such as - %"0.5 A + 0.5 B =>" by coincidence, when its magnitudes summed to - %exactly 1 --- and even then mishandled it below, since find() on such - %a reaction returns more than one index and only the first metabolite - %was ever used to build the boundary metabolite, silently dropping the - %rest. + %products. The coefficient magnitudes do not matter, so a scaled sink + %such as "2 A =>" and a multi-metabolite one-sided reaction such as + %"0.5 A + 0.5 B =>" are both boundary reactions. if any(col>0) && any(col<0) continue end @@ -53,7 +47,7 @@ continue end %A single metabolite names and inherits properties from that - %metabolite, as before. Several metabolites on the one populated side + %metabolite. Several metabolites on the one populated side %have no single metabolite to name the boundary after or copy %properties from, so it is named after the reaction instead and left %without those properties. diff --git a/src/manipulation/setExchangeBounds.m b/src/manipulation/setExchangeBounds.m index dddc6d5c..a118a175 100755 --- a/src/manipulation/setExchangeBounds.m +++ b/src/manipulation/setExchangeBounds.m @@ -200,10 +200,8 @@ % check for any metabolites that are exchanged in more than one reaction. % metInd holds one entry per (metabolite,reaction) pair, so a metabolite -% involved in more than one exchange reaction is a repeated *value* in -% metInd, not a position past numel(unique(metInd)) -- comparing those two -% index spaces directly (as a previous version of this check did) compares -% positions in the unique list against indices into metInd itself. +% involved in more than one exchange reaction is a value that occurs more +% than once in metInd. The occurrences are counted per unique value. [uniqueMetInd,~,ic] = unique(metInd); counts = accumarray(ic,1); multiMetInd = exchMetInd(uniqueMetInd(counts>1)); diff --git a/src/manipulation/simplifyModel.m b/src/manipulation/simplifyModel.m index 464f4011..2031936b 100755 --- a/src/manipulation/simplifyModel.m +++ b/src/manipulation/simplifyModel.m @@ -90,9 +90,10 @@ end if deleteDuplicates==true - %Delete all but the last occurrence of duplicate reactions. The - %reactions must have the same bounds, reversibility, and objective - %coefficient to be regarded as duplicate + %Merge duplicate reactions into the first occurrence and delete the + %others. Reactions with the same stoichiometry and reversibility are + %duplicates; the merged reaction keeps the widest bounds and the largest + %objective coefficient (see contractModel) [reducedModel, rxnsToDelete, ~]=contractModel(reducedModel); deletedReactions=[deletedReactions; rxnsToDelete]; end diff --git a/src/queries/getIndexes.m b/src/queries/getIndexes.m index 61840547..71fb1c59 100755 --- a/src/queries/getIndexes.m +++ b/src/queries/getIndexes.m @@ -25,7 +25,8 @@ % -------------------- % returnLogical : logical % sets whether to return a logical array or an array with the indexes -% (default false). +% (default false). It cannot be true for type 'metnames', which raises +% the error RAVEN:badInput. % % Returns % ------- @@ -61,6 +62,12 @@ case 'genes' searchIn=model.genes; case 'metnames' + if returnLogical==true + error('RAVEN:badInput', ['The combination of type ''metnames'' and ' ... + 'returnLogical true is not supported, because each name can ' ... + 'match several metabolites. Use type ''mets'' or ''metcomps'' ' ... + 'for a logical vector.']); + end searchIn=model.metNames; case 'metcomps' % If provided as metaboliteName[comp], then index search is quite diff --git a/src/queries/printFluxes.m b/src/queries/printFluxes.m index fde3fb9a..41654d17 100755 --- a/src/queries/printFluxes.m +++ b/src/queries/printFluxes.m @@ -46,8 +46,8 @@ function printFluxes(model, fluxes, varargin) % % Examples % -------- -% printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile, ... -% outputString, metaboliteList); +% printFluxes(model, fluxes, 'onlyExchange', false, 'cutOffFlux', 10^-6, ... +% 'outputFile', 'fluxes.txt', 'metaboliteList', {'glucose'}); p=parseRAVENargs(varargin, {'onlyExchange',true; 'cutOffFlux',10^-8; 'outputFile',[]; 'outputString',[]; 'metaboliteList',{}}); onlyExchange=p.onlyExchange; diff --git a/src/reconstruction/kegg/getModelFromKEGG.m b/src/reconstruction/kegg/getModelFromKEGG.m index 1b853173..f8712c21 100755 --- a/src/reconstruction/kegg/getModelFromKEGG.m +++ b/src/reconstruction/kegg/getModelFromKEGG.m @@ -24,7 +24,9 @@ % Returns % ------- % model : struct -% a model structure with all reactions and the metabolites used in them. +% a model structure with the reactions that the keep* options retain (by +% default all except those labelled as "general reaction") and the +% metabolites used in them. % KOModel : struct % a model structure representing KEGG Orthology ids and their associated % genes. KO ids are stored as reactions. diff --git a/tests/function_tests/tAnalysis.m b/tests/function_tests/tAnalysis.m index 9b49a086..e855d9e4 100644 --- a/tests/function_tests/tAnalysis.m +++ b/tests/function_tests/tAnalysis.m @@ -122,6 +122,24 @@ function findGeneDeletionsComputesGrRatio(testCase) testCase.verifyGreaterThan(max(grRatioMuts), 0); end + function findGeneDeletionsDoubleDeletionSizes(testCase) + % Three genes with solvable single deletions give three gene pairs. + m = struct(); + m.id = 'toy'; m.name = 'toy'; + m.rxns = {'EX_a'; 'R1'; 'R2'; 'R3'; 'BIO'; 'EX_p'}; m.rxnNames = m.rxns; + m.mets = {'a'; 'b'; 'p'}; m.metNames = m.mets; m.metComps = ones(3,1); + m.comps = {'c'}; m.compNames = {'c'}; + m.S = sparse([1 -1 -1 -1 0 0; 0 1 1 0 -1 0; 0 0 0 1 0 -1]); + m.lb = zeros(6,1); m.ub = [10; 1000; 1000; 1000; 1000; 1000]; + m.rev = zeros(6,1); m.c = [0; 0; 0; 0; 1; 0]; m.b = zeros(3,1); + m.genes = {'g1'; 'g2'; 'g3'}; m.grRules = {''; 'g1'; 'g2'; 'g3'; ''; ''}; + m.rxnGeneMat = sparse([0 0 0; 1 0 0; 0 1 0; 0 0 1; 0 0 0; 0 0 0]); + evalc('[genes, fluxes, ~, ~, grRatioMuts] = findGeneDeletions(m, ''dgd'');'); + testCase.verifySize(genes, [3 2]); + testCase.verifySize(fluxes, [numel(m.rxns) 3]); + testCase.verifySize(grRatioMuts, [1 3]); + end + function traceFluxPathReturnsCells(testCase) testCase.assumeSolver('solveLP'); sol = solveLP(testCase.model); diff --git a/tests/function_tests/tGapfilling.m b/tests/function_tests/tGapfilling.m index b4492fa6..04533ab8 100644 --- a/tests/function_tests/tGapfilling.m +++ b/tests/function_tests/tGapfilling.m @@ -26,6 +26,28 @@ function checkProductionReturnsIndices(testCase) testCase.verifyClass(notProduced, 'double'); end + function checkProductionReturnsExcretionReactionsOnly(testCase) + % Neither a nor c has a source, so b and d, and a and c, cannot be + % produced. The returned model holds an excretion reaction for + % each metabolite, with or without checkNeededForProduction. + m = struct(); + m.id = 'toy'; m.name = 'toy'; + m.rxns = {'R1'; 'R2'}; m.rxnNames = m.rxns; + m.mets = {'a'; 'b'; 'c'; 'd'}; m.metNames = m.mets; m.metComps = ones(4,1); + m.comps = {'c'}; m.compNames = {'c'}; + m.S = sparse([-1 0; 1 0; 0 -1; 0 1]); + m.lb = [0; 0]; m.ub = [1000; 1000]; m.rev = [0; 0]; m.c = [0; 0]; + m.b = zeros(4,1); + m.genes = {}; m.grRules = {''; ''}; m.rxnGeneMat = sparse(2,0); + + [~, ~, ~, ~, plain] = checkProduction(m, 'printDetails', false); + [notProduced, ~, ~, ~, checked] = checkProduction(m, ... + 'checkNeededForProduction', true, 'printDetails', false); + testCase.verifyEqual(notProduced, (1:4)'); + testCase.verifyEqual(full(checked.S), full(plain.S)); + testCase.verifyTrue(all(full(sum(checked.S(:, 3:end), 1)) < 0)); + end + function findLeakMetaboliteProduceReturnsSolution(testCase) evalc('[sol, metabolite] = findLeakMetabolite(testCase.model, ''produce'');'); testCase.verifyNotEmpty(sol); diff --git a/tests/function_tests/tGenomeData.m b/tests/function_tests/tGenomeData.m index 06855ece..d1dddff7 100644 --- a/tests/function_tests/tGenomeData.m +++ b/tests/function_tests/tGenomeData.m @@ -77,6 +77,25 @@ function renameModelGenesRenamesGenesAndGrRules(testCase) testCase.verifyEqual(size(m2.rxnGeneMat,2), numel(m2.genes)); end + function renameModelGenesReadsTableFromFilePath(testCase) + model.id = 'test'; + model.rxns = {'r1'; 'r2'}; + model.genes = {'LT1'; 'LT2'}; + model.grRules = {'LT1'; 'LT1 or LT2'}; + model.rxnGeneMat = sparse([1 0; 1 1]); + + T = table({'LT1'; 'LT2'}, {'AAA'; 'BBB'}, ... + 'VariableNames', {'locus_tag', 'gene_name'}); + tsvFile = [tempname '.tsv']; + writetable(T, tsvFile, 'FileType', 'text', 'Delimiter', '\t'); + testCase.addTeardown(@() delete(tsvFile)); + output = evalc('m2 = renameModelGenes(model, tsvFile, ''locus_tag'', ''gene_name'');'); + + testCase.verifySubstring(output, tsvFile); + testCase.verifyEqual(sort(m2.genes), {'AAA'; 'BBB'}); + testCase.verifyEqual(m2.grRules, {'AAA'; 'AAA or BBB'}); + end + function downloadGenomeDataDownloads(testCase) % Network test: download a small genome from NCBI Datasets. % Skipped when there is no internet connection. diff --git a/tests/function_tests/tQueries.m b/tests/function_tests/tQueries.m index 7e918050..7b146503 100644 --- a/tests/function_tests/tQueries.m +++ b/tests/function_tests/tQueries.m @@ -322,6 +322,19 @@ function getIndexesLogical(testCase) testCase.verifyEqual(find(lg), 1); end + function getIndexesMetnamesReturnLogicalErrors(testCase) + name = testCase.model.metNames(1); + testCase.verifyError( ... + @() getIndexes(testCase.model, name, 'metnames', 'returnLogical', true), ... + 'RAVEN:badInput'); + try + getIndexes(testCase.model, name, 'metnames', 'returnLogical', true); + catch e + testCase.verifySubstring(e.message, 'metnames'); + testCase.verifySubstring(e.message, 'returnLogical'); + end + end + function getIndexesLogicalMaskAllTrueReturnsIndices(testCase) % A logical all-true mask of length n must return numeric [1..n], % not be passed through as a logical array: this requires an