From 8bf11f4daa0698154ba34b6f6071d4459072e029 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:32:20 +0100 Subject: [PATCH 01/31] fillCarbonNum: warn and weight 1 for a model without metabolite formulas A model without metFormulas gave an error from the field access; every exchange now gets the carbon number 1 and a warning, as geckopy does. --- src/kcat_tuning/fillCarbonNum.m | 10 +++++- test/VerificationMatrix.tsv | 3 +- test/unit_tests/parityRound5G4Tests.m | 52 +++++++++++++++++++++++++++ 3 files changed, 63 insertions(+), 2 deletions(-) create mode 100644 test/unit_tests/parityRound5G4Tests.m diff --git a/src/kcat_tuning/fillCarbonNum.m b/src/kcat_tuning/fillCarbonNum.m index c976ea21..1e177bbe 100644 --- a/src/kcat_tuning/fillCarbonNum.m +++ b/src/kcat_tuning/fillCarbonNum.m @@ -2,7 +2,10 @@ % fillCarbonNum Add a carbon number field to the model. % % Adds a model.excarbon field with the carbon number of each exchange -% reaction, used to weight the RMSE during evotune kcat tuning. +% reaction, used to weight the RMSE during evotune kcat tuning. An exchange +% whose metabolite has no usable formula gets the carbon number 1. A model +% without a metFormulas field, or with only empty formulas, gives every +% exchange the carbon number 1 and a warning (GECKO:fillCarbonNum:noFormulas). % % Parameters % ---------- @@ -17,6 +20,11 @@ % ecModel with an added model.excarbon field. % get carbonnum for each exchange rxn to further calculation of error +if ~isfield(model,'metFormulas') || all(cellfun(@isempty,model.metFormulas)) + warning('GECKO:fillCarbonNum:noFormulas', ... + 'The model has no metabolite formulas; every exchange reaction gets the carbon number 1.') + model.metFormulas = repmat({''},numel(model.mets),1); +end model.excarbon = zeros(length(model.rxns),1); [EXrxn, EXrxnIdx] = getExchangeRxns(model); [CarbonNum,EXfors] = getcarbonnum(model,EXrxn); diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index bf792183..07346755 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -51,8 +51,9 @@ get_enzyme_data/getECstring.m Called from findECInDB - tested when that function download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 download_databases/fetchBrendaAsset.m tc0085 download_databases/loadDatabases.m Called from getECfromDatabase - tested when that function is tested. +kcat_tuning/evotune/bestParsimonious.m tc0600 kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases -kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test) +kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only), tc0601, tc0602 kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253, tc0254 kcat_tuning/sigmaFitter.m tc0020, tc0255 diff --git a/test/unit_tests/parityRound5G4Tests.m b/test/unit_tests/parityRound5G4Tests.m new file mode 100644 index 00000000..598a27b2 --- /dev/null +++ b/test/unit_tests/parityRound5G4Tests.m @@ -0,0 +1,52 @@ +%run this test case with the command +%results = runtests('parityRound5G4Tests.m') + +function tests = parityRound5G4Tests + tests = functiontests(localfunctions); +end + +function testBestParsimoniousSkipsNaNRmse_tc0600(testCase) + % A NaN rmse never sets the best rmse and is never chosen. + points = struct('threshold',{0,3,6}, 'nChanged',{4000,3000,1000}, ... + 'rmse',{NaN,0.905,1.4}); + point = bestParsimonious(points); + verifyEqual(testCase, point.nChanged, 3000) +end + +function testFillCarbonNumWithoutFormulasWarnsAndWeightsOne_tc0601(testCase) + model = exchangeToyModel(); + verifyWarning(testCase, @() fillCarbonNum(model), 'GECKO:fillCarbonNum:noFormulas') + model.metFormulas = {'';''}; + verifyWarning(testCase, @() fillCarbonNum(model), 'GECKO:fillCarbonNum:noFormulas') + warning('off','GECKO:fillCarbonNum:noFormulas') + cleanup = onCleanup(@() warning('on','GECKO:fillCarbonNum:noFormulas')); + withoutField = rmfield(model,'metFormulas'); + out = fillCarbonNum(withoutField); + verifyEqual(testCase, out.excarbon(strcmp(out.rxns,'EX_a')), 1) + verifyEqual(testCase, out.excarbon(strcmp(out.rxns,'EX_b')), 1) +end + +function testFillCarbonNumUsesFormulaCarbons_tc0602(testCase) + model = exchangeToyModel(); + model.metFormulas = {'C6H12O6';'O2'}; + out = fillCarbonNum(model); + verifyEqual(testCase, out.excarbon(strcmp(out.rxns,'EX_a')), 6) +end + +function model = exchangeToyModel() + model = struct(); + model.id = 'toy'; + model.rxns = {'EX_a';'EX_b'}; + model.rxnNames = {'exchange a';'exchange b'}; + model.mets = {'a';'b'}; + model.metNames = {'a';'b'}; + model.S = sparse([-1 0;0 -1]); + model.lb = [-1000;-1000]; + model.ub = [1000;1000]; + model.rev = [1;1]; + model.c = [0;0]; + model.b = [0;0]; + model.comps = {'c'}; + model.compNames = {'cytosol'}; + model.metComps = [1;1]; +end From c39063f8310f3ff4e0df3e6f87f12655e8f6a337 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:35:18 +0100 Subject: [PATCH 02/31] ecStatsFromBrenda gives the median and the maximum per EC and keeps the EC case --- src/gather_kcats/ecStatsFromBrenda.m | 69 +++++++++++++++++------- test/unit_tests/brendaReviewStatsTests.m | 53 ++++++++++++++++++ 2 files changed, 102 insertions(+), 20 deletions(-) create mode 100644 test/unit_tests/brendaReviewStatsTests.m diff --git a/src/gather_kcats/ecStatsFromBrenda.m b/src/gather_kcats/ecStatsFromBrenda.m index 45775ebf..e14b291d 100644 --- a/src/gather_kcats/ecStatsFromBrenda.m +++ b/src/gather_kcats/ecStatsFromBrenda.m @@ -1,41 +1,61 @@ -function ecStats = ecStatsFromBrenda(KCATcell,SAcell) +function ecStats = ecStatsFromBrenda(KCATcell,SAcell,KCATmedianCell) % ecStatsFromBrenda Summarise loadBRENDAdata's output per EC number, for % reviewKcatAssignment's ecStats argument. % -% TODO: verify against a real BRENDA pull before relying on this. Ported -% from geckopy's ec_stats_from_brenda without being able to run it -% against real BRENDA data in this environment. In particular, -% loadBRENDAdata's KCATcell only carries the BRENDA *max* aggregate per -% (EC, substrate, organism) triple (column 4 of kcat.tsv) -- the median -% column is read internally but discarded before KCATcell is returned -- -% so kcatMedian below is always NaN, and reviewKcatAssignment's -% "ec-maximum" check (which needs a finite kcatMedian) will never fire -% until loadBRENDAdata is extended to also expose it. +% loadBRENDAdata returns one aggregate of kcat.tsv per call, so the maximum +% view and the median view come from two calls: +% +% [KCATcell, SAcell] = loadBRENDAdata('kcatAggregation', 'max'); +% KCATmedianCell = loadBRENDAdata('kcatAggregation', 'median'); +% ecStats = ecStatsFromBrenda(KCATcell, SAcell, KCATmedianCell); +% +% geckopy's ec_stats_from_brenda gives the same statistics. The EC numbers +% keep the case they have in the BRENDA files. % % Parameters % ---------- % KCATcell : cell -% as returned by loadBRENDAdata's first output: {ecCode, substrate, -% organism, kcatMax}. +% first output of loadBRENDAdata with kcatAggregation 'max': {ecCode, +% substrate, organism, kcatMax}. % SAcell : cell -% as returned by loadBRENDAdata's second output: {ecCode, organism, -% kcatFromSA, mw}. +% second output of loadBRENDAdata with kcatAggregation 'max': {ecCode, +% organism, kcatFromSA, mw}. +% KCATmedianCell : cell +% first output of loadBRENDAdata with kcatAggregation 'median' +% (optional). Without it kcatMedian is NaN, so reviewKcatAssignment +% cannot flag 'ec-maximum', and a warning is given. % % Returns % ------- % ecStats : containers.Map -% EC code -> struct(nKcat,kcatMax,kcatMedian,values,nSa), keyed by -% the same (uppercased) EC strings BRENDA uses. +% EC code -> struct(nKcat,kcatMax,kcatMedian,values,nSa). nKcat and +% kcatMax are the number and the maximum of the kcat rows of the +% maximum view, kcatMedian is the median over the rows of the median +% view, values holds the kcat values of both views rounded to six +% decimals, and nSa is the number of specific-activity rows. % % See also % -------- % reviewKcatAssignment, loadBRENDAdata -ecCol = upper(KCATcell{1}); +ecCol = KCATcell{1}; maxCol = KCATcell{4}; +hasMedian = nargin >= 3 && ~isempty(KCATmedianCell) && ~isempty(KCATmedianCell{1}); +if nargin < 3 + warning('ecStatsFromBrenda:noMedianView', ... + ['No median view of the kcat rows was given, so kcatMedian is NaN ' ... + 'and reviewKcatAssignment cannot flag ec-maximum.']) +end +medEcCol = {}; +medCol = []; +if hasMedian + medEcCol = KCATmedianCell{1}; + medCol = KCATmedianCell{4}; +end + saEcCol = {}; if ~isempty(SAcell) && numel(SAcell) >= 1 && ~isempty(SAcell{1}) - saEcCol = upper(SAcell{1}); + saEcCol = SAcell{1}; end ecList = unique([ecCol(:); saEcCol(:)]); @@ -46,12 +66,21 @@ nKcat = sum(rows); if nKcat > 0 kcatMax = max(maxCol(rows)); - values = unique(round(maxCol(rows),6)); + values = round(maxCol(rows),6); else kcatMax = NaN; values = []; end + kcatMedian = NaN; + if hasMedian + medRows = strcmp(medEcCol,ec); + if any(medRows) + kcatMedian = median(medCol(medRows),'omitnan'); + values = [values(:); round(medCol(medRows),6)]; + end + end + values = unique(values(:)); nSa = sum(strcmp(saEcCol,ec)); - ecStats(ec) = struct('nKcat',nKcat,'kcatMax',kcatMax,'kcatMedian',NaN,'values',values,'nSa',nSa); + ecStats(ec) = struct('nKcat',nKcat,'kcatMax',kcatMax,'kcatMedian',kcatMedian,'values',values,'nSa',nSa); end end diff --git a/test/unit_tests/brendaReviewStatsTests.m b/test/unit_tests/brendaReviewStatsTests.m new file mode 100644 index 00000000..83b98203 --- /dev/null +++ b/test/unit_tests/brendaReviewStatsTests.m @@ -0,0 +1,53 @@ +%run this test case with the command +%results = runtests('brendaReviewStatsTests.m') + +function tests = brendaReviewStatsTests + tests = functiontests(localfunctions); +end + +function testEcStatsFromBrendaMedianMaximumAndCase_tc0400(testCase) + % Per EC number: the rows and maximum of the max view, the median of the + % median view, the values of both views, and the specific-activity rows. + % The EC keys keep their case. + kcatMax = {{'1.1.1.1';'1.1.1.1';'3.5.1.n3';'3.5.1.N5'}, {'a';'b';'c';'d'}, {'o';'o';'o';'o'}, [10;4;40;7]}; + kcatMed = {{'1.1.1.1';'1.1.1.1';'1.1.1.1';'3.5.1.n3'}, {'a';'b';'e';'c'}, {'o';'o';'o';'o'}, [2;3;6;4]}; + sa = {{'4.1.1.1';'1.1.1.1'}, {'o';'o'}, [1;2], [1;1]}; + stats = ecStatsFromBrenda(kcatMax, sa, kcatMed); + verifyEqual(testCase, sort(stats.keys()), {'1.1.1.1','3.5.1.N5','3.5.1.n3','4.1.1.1'}) + s = stats('1.1.1.1'); + verifyEqual(testCase, [s.nKcat s.kcatMax s.kcatMedian s.nSa], [2 10 3 1]) + verifyEqual(testCase, s.values(:)', [2 3 4 6 10]) + s = stats('3.5.1.n3'); + verifyEqual(testCase, [s.nKcat s.kcatMax s.kcatMedian s.nSa], [1 40 4 0]) + verifyEqual(testCase, s.values(:)', [4 40]) + s = stats('3.5.1.N5'); + verifyEqual(testCase, [s.nKcat s.kcatMax s.nSa], [1 7 0]) + verifyTrue(testCase, isnan(s.kcatMedian)) + s = stats('4.1.1.1'); + verifyEqual(testCase, [s.nKcat s.nSa], [0 1]) + verifyTrue(testCase, isnan(s.kcatMax) && isnan(s.kcatMedian)) + verifyEmpty(testCase, s.values) +end + +function testEcStatsFromBrendaWithoutMedianViewWarns_tc0401(testCase) + kcatMax = {{'1.1.1.1';'1.1.1.1'}, {'a';'b'}, {'o';'o'}, [10;4]}; + verifyWarning(testCase, @() ecStatsFromBrenda(kcatMax, {[],[],[],[]}), 'ecStatsFromBrenda:noMedianView') + warning('off', 'ecStatsFromBrenda:noMedianView') + cleanup = onCleanup(@() warning('on', 'ecStatsFromBrenda:noMedianView')); %#ok + stats = ecStatsFromBrenda(kcatMax, {[],[],[],[]}); + s = stats('1.1.1.1'); + verifyTrue(testCase, isnan(s.kcatMedian)) + verifyEqual(testCase, [s.nKcat s.kcatMax], [2 10]) + verifyEqual(testCase, s.values(:)', [4 10]) +end + +function testReviewFlagsEcMaximumWithLoaderStatistics_tc0402(testCase) + % A kcat equal to the maximum of its EC, ten times the median, is flagged + % 'ec-maximum', also for an EC with a lowercase letter. + kcatMax = {{'1.1.1.1';'3.5.1.n3'}, {'a';'c'}, {'o';'o'}, [10;40]}; + kcatMed = {{'1.1.1.1';'3.5.1.n3'}, {'a';'c'}, {'o';'o'}, [2;4]}; + stats = ecStatsFromBrenda(kcatMax, {[],[],[],[]}, kcatMed); + f = reviewKcatAssignment({'R_A';'R_F'}, [10;40], {'brenda';'brenda'}, {'1.1.1.1';'3.5.1.n3'}, 'ecStats', stats); + verifyEqual(testCase, f.rxnId', {'R_A','R_F'}) + verifyTrue(testCase, all(contains(f.checks, 'ec-maximum'))) +end From a71a026be5507afff4529cece98607c23f7475e7 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:35:18 +0100 Subject: [PATCH 03/31] loadBrendaMeasurements reads measurements.tsv.xz --- src/gather_kcats/loadBrendaMeasurements.m | 28 ++++++++++++-- test/unit_tests/brendaMeasurementsXzTests.m | 42 +++++++++++++++++++++ 2 files changed, 66 insertions(+), 4 deletions(-) create mode 100644 test/unit_tests/brendaMeasurementsXzTests.m diff --git a/src/gather_kcats/loadBrendaMeasurements.m b/src/gather_kcats/loadBrendaMeasurements.m index 52d66913..3023f4c7 100644 --- a/src/gather_kcats/loadBrendaMeasurements.m +++ b/src/gather_kcats/loadBrendaMeasurements.m @@ -1,8 +1,10 @@ function measurements = loadBrendaMeasurements(varargin) % loadBrendaMeasurements Load the BRENDA measurement table. % -% Reads measurements.tsv from the BRENDA database folder of the model -% adapter (next to kcat.tsv, sa.tsv and mw.tsv): one row per wild-type +% Reads measurements.tsv.xz, or measurements.tsv when there is no .xz file, +% from the BRENDA database folder of the model adapter (next to kcat.tsv, +% sa.tsv and mw.tsv). A .xz file is decompressed with the xz command, which +% has to be on the system path. The table has one row per wild-type % BRENDA kcat with an EC number, after the outlier screen, with the % substrate, organism, UniProt accession and sequence id of each % measurement. geckopy reads the same file (load_brenda_measurements). @@ -12,7 +14,8 @@ % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % folder : char -% folder with measurements.tsv, instead of the adapter's BRENDA folder. +% folder with measurements.tsv.xz or measurements.tsv, instead of the +% adapter's BRENDA folder. % % Returns % ------- @@ -43,8 +46,19 @@ end folder = modelAdapter.getBrendaDBFolder(); end +xzFile = fullfile(folder, 'measurements.tsv.xz'); file = fullfile(folder, 'measurements.tsv'); -if ~isfile(file) +tempFile = ''; +if isfile(xzFile) + tempFile = [tempname() '.tsv']; + cleanup = onCleanup(@() deleteIfExists(tempFile)); %#ok + [status, output] = system(sprintf('xz -dc "%s" > "%s"', xzFile, tempFile)); + if status ~= 0 + error('loadBrendaMeasurements:xz', ['Could not decompress %s with the xz command ' ... + '(install xz or decompress the file to measurements.tsv): %s'], xzFile, strtrim(output)) + end + file = tempFile; +elseif ~isfile(file) error('loadBrendaMeasurements:missing', ['No BRENDA measurement table (measurements.tsv) in %s. ' ... 'It comes with the BRENDA files of GECKO (databases folder).'], folder) end @@ -75,6 +89,12 @@ end end +function deleteIfExists(file) +if isfile(file) + delete(file) +end +end + function checkOptions(args, names, functionName) % Options are name-value pairs with known names, so a misspelled option is % an error instead of being taken as a positional argument. diff --git a/test/unit_tests/brendaMeasurementsXzTests.m b/test/unit_tests/brendaMeasurementsXzTests.m new file mode 100644 index 00000000..7a9f361d --- /dev/null +++ b/test/unit_tests/brendaMeasurementsXzTests.m @@ -0,0 +1,42 @@ +%run this test case with the command +%results = runtests('brendaMeasurementsXzTests.m') + +function tests = brendaMeasurementsXzTests + tests = functiontests(localfunctions); +end + +function testLoadBrendaMeasurementsReadsXz_tc0403(testCase) + [status, ~] = system('xz --version'); + assumeEqual(testCase, status, 0) + folder = tempname(); + mkdir(folder); + cleanup = onCleanup(@() rmdir(folder, 's')); %#ok + writeMeasurements(fullfile(folder, 'measurements.tsv'), {'1.1.1.1', 5}) + % only the plain file + m = loadBrendaMeasurements('folder', folder); + verifyEqual(testCase, m.kcat, 5) + % only the compressed file + verifyEqual(testCase, system(sprintf('xz -k "%s"', fullfile(folder, 'measurements.tsv'))), 0) + delete(fullfile(folder, 'measurements.tsv')) + m = loadBrendaMeasurements('folder', folder); + verifyEqual(testCase, m.kcat, 5) + verifyEqual(testCase, m.release, 'BRENDA release 2026.1') + % both files: the compressed one is read + writeMeasurements(fullfile(folder, 'measurements.tsv'), {'1.1.1.1', 7}) + m = loadBrendaMeasurements('folder', folder); + verifyEqual(testCase, m.kcat, 5) + % neither + delete(fullfile(folder, 'measurements.tsv')) + delete(fullfile(folder, 'measurements.tsv.xz')) + verifyError(testCase, @() loadBrendaMeasurements('folder', folder), 'loadBrendaMeasurements:missing') +end + +function writeMeasurements(file, rows) + fid = fopen(file, 'w'); + fprintf(fid, '# BRENDA release 2026.1\n'); + fprintf(fid, 'ec\tsubstrate\tsubstrate_inchikey\torganism\ttaxid\tuniprot\tseq_id\tkcat\ttemperature\tph\tpmid\n'); + for i = 1:size(rows, 1) + fprintf(fid, '%s\tS\tK\tO\t1\tP1\t1\t%g\t25\t7\t1\n', rows{i, 1}, rows{i, 2}); + end + fclose(fid); +end From cc6b106e4a344accb9f3cc848e8dc43b1a2ffd2f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:35:18 +0100 Subject: [PATCH 04/31] getECfromDatabase warns on multiplicities with action display and searches the selected reactions only --- src/get_enzyme_data/getECfromDatabase.m | 22 ++++++--- test/unit_tests/ecFromDatabaseConflictTests.m | 45 +++++++++++++++++++ 2 files changed, 61 insertions(+), 6 deletions(-) create mode 100644 test/unit_tests/ecFromDatabaseConflictTests.m diff --git a/src/get_enzyme_data/getECfromDatabase.m b/src/get_enzyme_data/getECfromDatabase.m index 3fd35fc7..16e88f3f 100644 --- a/src/get_enzyme_data/getECfromDatabase.m +++ b/src/get_enzyme_data/getECfromDatabase.m @@ -16,15 +16,16 @@ % of length model.ec.rxns that specifies which model.ec.eccodes entries % should be queried. Existing values in model.ec.eccodes will be wiped. % Entries that are indicated by false will be kept and not modified by -% this function (by default all model.ec.eccodes entries are populated -% by this function). +% this function, and are not searched for multiplicities (by default all +% model.ec.eccodes entries are populated by this function). % action : char % response action if multiple proteins with different EC numbers are % found for a given gene in a metabolic reaction (default 'display'): % -% - 'display' displays all found multiplicities. -% - 'ignore' ignore multiplicities and use the protein with the lowest -% index in the database. +% - 'display' gives a warning that lists all found multiplicities and +% uses the protein with the lowest index in the database. +% - 'ignore' ignores multiplicities and uses the protein with the +% lowest index in the database. % - 'add' adds all the multiple proteins as isozymes for the given % reaction. % modelAdapter : ModelAdapter @@ -96,9 +97,17 @@ conflicts = cell(1,4); rxnEnzMat = logical(rxnEnzMat); +if isempty(ecRxns) + selected = true(n,1); +else + selected = logical(ecRxns(:)); +end PB = progressReport(n, 'Assigning EC numbers from database'); for i = 1:n PB.count; + if ~selected(i) + continue + end gns = genes(rxnEnzMat(i,:).'); if ~isempty(gns) %Find match in Uniprot: @@ -176,10 +185,11 @@ function displayErrorMessage(conflicts,uniprot,kegg) STR = [STR, '\nIf a wrongly annotated case was found then call the ']; STR = [STR, 'getECfromDatabase.m function again with the option action']; STR = [STR, '= ignore\n\n']; +STR = [STR, 'The protein with the lowest index in the database was used for each of them.\n\n']; STR = [STR, 'If the conflicting proteins are desired to be kept as isozymes']; STR = [STR, ' then call the getECfromDatabase.m function']; STR = [STR, ' again with the option action = add\n']; -error(sprintf(STR)) +warning('getECfromDatabase:conflicts', '%s', sprintf(STR)) end function [geneIndex,geneHashMap]=hashGeneToProt(proteinDB) diff --git a/test/unit_tests/ecFromDatabaseConflictTests.m b/test/unit_tests/ecFromDatabaseConflictTests.m new file mode 100644 index 00000000..a8cc6027 --- /dev/null +++ b/test/unit_tests/ecFromDatabaseConflictTests.m @@ -0,0 +1,45 @@ +%run this test case with the command +%results = runtests('ecFromDatabaseConflictTests.m') + +function tests = ecFromDatabaseConflictTests + tests = functiontests(localfunctions); +end + +function testGetECfromDatabaseDisplayWarnsAndUsesFirst_tc0404(testCase) + [ecModel, adapter] = conflictModel(testCase); + verifyWarning(testCase, @() getECfromDatabase(ecModel, [], 'display', adapter), 'getECfromDatabase:conflicts') + ecModel2 = getECfromDatabase(ecModel, [], 'ignore', adapter); + verifyEqual(testCase, ecModel2.ec.eccodes, {'1.1.1.1';'1.1.1.1';'1.1.1.1';'1.1.1.1';'1.1.2.1';'1.1.1.3'}) + warning('off', 'getECfromDatabase:conflicts') + cleanup = onCleanup(@() warning('on', 'getECfromDatabase:conflicts')); %#ok + displayed = getECfromDatabase(ecModel, [], 'display', adapter); + verifyEqual(testCase, displayed.ec.eccodes, ecModel2.ec.eccodes) +end + +function testGetECfromDatabaseIgnoresConflictOutsideSelection_tc0405(testCase) + [ecModel, adapter] = conflictModel(testCase); + selected = strcmp(ecModel.ec.rxns, 'R2_EXP_1'); + lastwarn('') + ecModel2 = getECfromDatabase(ecModel, selected, 'display', adapter); + verifyEmpty(testCase, lastwarn()) + verifyEqual(testCase, ecModel2.ec.eccodes(selected), {'1.1.1.1'}) + verifyEqual(testCase, ecModel2.ec.eccodes(~selected), repmat({''}, nnz(~selected), 1)) +end + +function [ecModel, adapter] = conflictModel(testCase) + % ecTestGEM with a second UniProt entry for gene G4 and another EC number. + geckoPath = findGECKOroot; + source = fullfile(geckoPath,'test','unit_tests','ecTestGEM'); + adapter = ModelAdapterManager.getAdapter(fullfile(source, 'TestGEMAdapter.m')); + model = getGeckoTestModel(); + ecModel = makeEcModel(model, false, adapter); + folder = tempname(); + mkdir(fullfile(folder, 'data')); + testCase.addTeardown(@() rmdir(folder, 's')) + copyfile(fullfile(source, 'data', 'uniprot.tsv'), fullfile(folder, 'data', 'uniprot.tsv')) + copyfile(fullfile(source, 'data', 'kegg.tsv'), fullfile(folder, 'data', 'kegg.tsv')) + fid = fopen(fullfile(folder, 'data', 'uniprot.tsv'), 'a'); + fprintf(fid, '\nP4b\tG4\t1.1.2.2\t40001\tMDFM\n'); + fclose(fid); + adapter.params.path = folder; +end From 3028b24839a4ee89bb3f95c2c4ae7a24046f850f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:35:18 +0100 Subject: [PATCH 05/31] Record the new tests in the verification matrix --- test/VerificationMatrix.tsv | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index bf792183..2ea31f48 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -26,7 +26,7 @@ gather_kcats/collectKcats.m tc0079, tc0080, tc0081, tc0083, tc0087, tc0088, tc00 gather_kcats/assignKcats.m tc0079, tc0083 gather_kcats/updateKcatSources.m tc0082, tc0087 gather_kcats/kcatAssignmentConfig.m tc0079, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090 -gather_kcats/loadBrendaMeasurements.m tc0085 +gather_kcats/loadBrendaMeasurements.m tc0085, tc0403 gather_kcats/ecTreeEstimate.m tc0084 gather_kcats/kcatsFromECTree.m tc0086, tc0087 gather_kcats/kcatsFromHomologs.m tc0095, tc0096 @@ -40,12 +40,12 @@ flux_analysis/ecFVA.m tc0214 flux_analysis/ecFSEOF.m tc0100, tc0236 utilities/upgradeGECKOproject.m tc0094, tc0099 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 -gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases -gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment +gather_kcats/reviewKcatAssignment.m tc0402 (the ec-maximum check); the other checks have only been tested by running it manually +gather_kcats/ecStatsFromBrenda.m tc0400, tc0401 - not verified against a real BRENDA pull gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243, tc0244, tc0245, tc0246 -get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM +get_enzyme_data/getECfromDatabase.m tc0007, tc0404, tc0405 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested; also tc0246 download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 From c1941dff7ec14c0a5bce37a0c483c2779cdad9c0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 06/31] addNewRxnsToEC adds reactions whose genes already have enzymes --- src/change_model/addNewRxnsToEC.m | 35 ++++++++++++++++++------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/src/change_model/addNewRxnsToEC.m b/src/change_model/addNewRxnsToEC.m index 2ba53466..9648356c 100644 --- a/src/change_model/addNewRxnsToEC.m +++ b/src/change_model/addNewRxnsToEC.m @@ -50,7 +50,9 @@ % Note that the reversibility defined in the equation will be used to % split to construct irreversible reactions (add _REV) and the rules % defined will be used to expand the model (add _EXP_n). -% 2. After adding the new reactions, setKcatForReactions or +% 2. When all the genes of the new reactions already have enzymes, no +% newEnzymes are needed. +% 3. After adding the new reactions, setKcatForReactions or % applyCustomKcats should be run to add the kcat to the reactions, and % subsequently applyKcatConstraints. % @@ -78,8 +80,9 @@ end params = modelAdapter.getParameters(); +enzAdded = {}; if isempty(newEnzymes) - newEnzymes = []; + newEnzymes = struct('enzymes', {{}}, 'genes', {{}}, 'mw', []); end % Validate if GECKO version 3 @@ -91,7 +94,7 @@ error('addNewRxnsToEC only works on non-light ecModels') else % Check if the new proteins are not already present in the model - if ~isempty(newEnzymes) + if ~isempty(newEnzymes.enzymes) [toRemove,~] = ismember(newEnzymes.enzymes, model.ec.enzymes); % Remove from the list those enzymes that already present @@ -213,19 +216,21 @@ model = addRxns(model, newRxns, 3, '', true, true); % 4. Add protein usage reactions. - usageRxns.rxns = strcat('usage_',proteinMets.mets); - usageRxns.rxnNames = usageRxns.rxns; - usageRxns.mets = cell(numel(usageRxns.rxns),1); - usageRxns.stoichCoeffs = cell(numel(usageRxns.rxns),1); - for i=1:numel(usageRxns.mets) - usageRxns.mets{i} = {proteinMets.mets{i}, 'prot_pool'}; - usageRxns.stoichCoeffs{i} = [1,-1]; + if ~isempty(newEnzymes.enzymes) + usageRxns.rxns = strcat('usage_',proteinMets.mets); + usageRxns.rxnNames = usageRxns.rxns; + usageRxns.mets = cell(numel(usageRxns.rxns),1); + usageRxns.stoichCoeffs = cell(numel(usageRxns.rxns),1); + for i=1:numel(usageRxns.mets) + usageRxns.mets{i} = {proteinMets.mets{i}, 'prot_pool'}; + usageRxns.stoichCoeffs{i} = [1,-1]; + end + usageRxns.lb = zeros(numel(usageRxns.rxns),1); + usageRxns.ub = zeros(numel(usageRxns.rxns),1) + 1000; + usageRxns.rev = ones(numel(usageRxns.rxns),1); + usageRxns.grRules = newEnzymes.genes; + model = addRxns(model,usageRxns); end - usageRxns.lb = zeros(numel(usageRxns.rxns),1); - usageRxns.ub = zeros(numel(usageRxns.rxns),1) + 1000; - usageRxns.rev = ones(numel(usageRxns.rxns),1); - usageRxns.grRules = newEnzymes.genes; - model = addRxns(model,usageRxns); % 5. Update .ec field with enzyme data model.ec.enzymes = [model.ec.enzymes; newEnzymes.enzymes(:)]; From 7ee7af795aef78e348c33e7e6f1b61dd76995d2e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 07/31] applyComplexData counts a subunit without stoichiometry as one subunit --- src/change_model/applyComplexData.m | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/change_model/applyComplexData.m b/src/change_model/applyComplexData.m index 68bcc976..2cc17579 100644 --- a/src/change_model/applyComplexData.m +++ b/src/change_model/applyComplexData.m @@ -38,7 +38,9 @@ % Notes % ----- % The subunit stoichiometry is assigned per protein ID: the proteins of a -% complex need not be listed in the order of model.ec.enzymes. +% complex need not be listed in the order of model.ec.enzymes. A subunit with +% stoichiometry 0 has no stated stoichiometry; it counts as a subunit with +% stoichiometry 1. % % foundComplex and proposedComplex share the following columns: rxn % (reaction identifier), complexID and name (Complex Portal identifiers), @@ -113,9 +115,8 @@ protMatIdx(~isOld) = lastProt+1 : lastProt+numel(newProts); lastProt = lastProt + numel(newProts); % Add subunit numbers - if all(complexData(i).stochiometry == 0) - complexData(i).stochiometry(:) = 1; - end + % A subunit without a stated stoichiometry (0) counts as one subunit + complexData(i).stochiometry(complexData(i).stochiometry == 0) = 1; complexMatrix(i,protMatIdx) = complexData(i).stochiometry; end From 76238c7d6251fe566feb3e9476b8bac0b49a333a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 08/31] calculateFfactor and relaxProteomicsGreedy read an empty cell array as no enzymes --- src/limit_proteins/calculateFfactor.m | 5 +++-- src/limit_proteins/relaxProteomicsGreedy.m | 6 ++++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/limit_proteins/calculateFfactor.m b/src/limit_proteins/calculateFfactor.m index 16796059..7cac14fc 100644 --- a/src/limit_proteins/calculateFfactor.m +++ b/src/limit_proteins/calculateFfactor.m @@ -16,7 +16,8 @@ % attempts to load data/paxDB.tsv). With several samples, the abundance % of a protein is averaged over the samples in which it is measured. % enzymes : cell -% list of enzymes (default model.ec.enzymes). +% list of enzymes (default model.ec.enzymes). An empty cell array {} +% means no enzymes (f is 0); [] means the default. % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % @@ -51,7 +52,7 @@ end params = modelAdapter.getParameters(); -if isempty(enzymes) +if isnumeric(enzymes) && isempty(enzymes) enzymes = model.ec.enzymes; end diff --git a/src/limit_proteins/relaxProteomicsGreedy.m b/src/limit_proteins/relaxProteomicsGreedy.m index 91a3df49..60171c3c 100644 --- a/src/limit_proteins/relaxProteomicsGreedy.m +++ b/src/limit_proteins/relaxProteomicsGreedy.m @@ -19,7 +19,9 @@ % target objective value. Required. % enzymeSet : cell % subset of uniprot ids eligible for relaxation (default: every -% currently proteomics-constrained enzyme). +% currently proteomics-constrained enzyme). An empty cell array {} means +% no enzyme is eligible (converged is false unless the model already +% reaches minimalGrowth); [] means the default. % maxIterations : double % safety cap; errors if exceeded while eligible enzymes still remained % (default 100). @@ -91,7 +93,7 @@ error('Cannot find prot_ metabolites for all enzymes. This is done by makeEcModel.') end -if isempty(enzymeSet) +if isnumeric(enzymeSet) && isempty(enzymeSet) pool = model.ec.enzymes; else pool = intersect(enzymeSet(:), model.ec.enzymes, 'stable'); From 77e7dd6298fe1c57867151f8f4a9e807e3361c58 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 09/31] constrainEnzConcs and fillEnzConcs raise for a negative concentration --- src/limit_proteins/constrainEnzConcs.m | 8 +++++++- src/limit_proteins/fillEnzConcs.m | 8 ++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/limit_proteins/constrainEnzConcs.m b/src/limit_proteins/constrainEnzConcs.m index 5029b400..75b90385 100644 --- a/src/limit_proteins/constrainEnzConcs.m +++ b/src/limit_proteins/constrainEnzConcs.m @@ -25,7 +25,8 @@ % % Notes % ----- -% To populate model.ec.concs you should run fillEnzConcs. +% A negative concentration in model.ec.concs is contradictory data and +% raises the error constrainEnzConcs:negativeConcentration. To populate model.ec.concs you should run fillEnzConcs. % % Examples % -------- @@ -63,6 +64,11 @@ protCons = []; else protCons = ~isnan(model.ec.concs); + if any(model.ec.concs(protCons) < 0) + error('constrainEnzConcs:negativeConcentration', ... + 'model.ec.concs has a negative concentration for %s.', ... + strjoin(model.ec.enzymes(model.ec.concs < 0), ', ')) + end end %Set all reactions to draw from prot_pool (forward: prot_pool -> prot_) diff --git a/src/limit_proteins/fillEnzConcs.m b/src/limit_proteins/fillEnzConcs.m index 501567ed..6353d012 100644 --- a/src/limit_proteins/fillEnzConcs.m +++ b/src/limit_proteins/fillEnzConcs.m @@ -32,6 +32,9 @@ % % Notes % ----- +% A negative abundance of a protein that is an enzyme of the model is +% contradictory data and raises the error fillEnzConcs:negativeAbundance. +% % To also constrain the model with the content of model.ec.concs, you should % run constrainEnzConcs. % @@ -57,5 +60,10 @@ model.ec.concs=nan(numel(model.ec.enzymes),1); [a,b] = ismember(uniprotIDs, model.ec.enzymes); +if any(abundances(a) < 0) + error('fillEnzConcs:negativeAbundance', ... + 'protData has a negative abundance for %s.', ... + strjoin(uniprotIDs(a & abundances < 0), ', ')) +end model.ec.concs(b(a)) = abundances(a); end From 6f773b5d31d8192dfac1d2971a97c88e82d9228a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 10/31] constrainFluxData raises for an invalid maxMinGrowth, a missing reaction and an unknown condition --- src/limit_proteins/constrainFluxData.m | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/src/limit_proteins/constrainFluxData.m b/src/limit_proteins/constrainFluxData.m index fc10d369..fa998c05 100644 --- a/src/limit_proteins/constrainFluxData.m +++ b/src/limit_proteins/constrainFluxData.m @@ -27,7 +27,8 @@ % 'max' if the provided growth rate should be set as maximum growth rate % (= upper bound), or 'min' if it should be set as minimum growth rate % (= lower bound). The latter option is suitable if minimization of -% prot_pool_exchange is used as objective function (default 'max'). +% prot_pool_exchange is used as objective function (default 'max'). Any +% other value raises constrainFluxData:invalidMaxMinGrowth. % looseStrictFlux : char or double % how strictly constrained the exchange fluxes should be (default % 'loose'): @@ -49,6 +50,10 @@ % % Notes % ----- +% An exchange reaction in fluxData.exchRxnIDs that is not in the model +% raises constrainFluxData:missingReaction, and a condition that is not in +% fluxData.conds raises constrainFluxData:unknownCondition. +% % If a provided constraint is either -1000 or 1000, then the function will % update the reaction lower and upper bound to either allow uptake or % excretion, irrespective of what option is given as the looseStrictFlux @@ -96,12 +101,24 @@ fluxData = loadFluxData(fullfile(params.path,'data','fluxData.tsv'),modelAdapter); end +if ~(ischar(maxMinGrowth) && any(strcmp(maxMinGrowth, {'max','min'}))) + error('constrainFluxData:invalidMaxMinGrowth', ... + 'maxMinGrowth must be ''max'' or ''min''.') +end + +missingRxns = fluxData.exchRxnIDs(~ismember(fluxData.exchRxnIDs, model.rxns)); +if ~isempty(missingRxns) + error('constrainFluxData:missingReaction', ... + 'The exchange reactions %s of fluxData are not in the model.', strjoin(missingRxns, ', ')) +end + if isempty(condition) condition = 1; elseif ~isnumeric(condition) idx = find(strcmp(fluxData.conds,condition)); - if isempty(condition) - error(['Condition ' condition ' cannot be found in fluxData']) + if isempty(idx) + error('constrainFluxData:unknownCondition', ... + 'Condition %s cannot be found in fluxData.', condition) else condition = idx; end From e425e5ea1a24eb3a0948f3e47a17c70e557e1b42 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 11/31] loadFluxData raises for a missing file, a ragged row and no exchange columns, and reads nested parentheses in headers --- src/io/loadFluxData.m | 63 +++++++++++++++++++---- src/kcat_tuning/evotune/loadEvotuneData.m | 17 ++++-- 2 files changed, 64 insertions(+), 16 deletions(-) diff --git a/src/io/loadFluxData.m b/src/io/loadFluxData.m index 49ec62cb..bd14a404 100644 --- a/src/io/loadFluxData.m +++ b/src/io/loadFluxData.m @@ -19,6 +19,13 @@ % % Notes % ----- +% The function raises an error when the file does not exist +% (loadFluxData:fileNotFound), when a line has a different number of cells +% than the header (loadFluxData:raggedRow) or when the file has no exchange +% flux columns (loadFluxData:noExchangeColumns). An exchange column header has +% the form ' ()'; the reaction id may contain +% parentheses. +% % The fluxData structure has the following fields: % % - conds : sampling condition. @@ -59,21 +66,23 @@ end %Load total protein content and flux data -try - fID = fopen(fluxDataFile); - formatStr = '%s'; - data = textscan(fID,formatStr,'Delimiter','\n'); - fclose(fID); -catch - warning('Failed to load %s.',fluxDataFile); - fluxData = []; - return +if ~isfile(fluxDataFile) + error('loadFluxData:fileNotFound', 'The flux data file %s does not exist.', fluxDataFile); end +fID = fopen(fluxDataFile); +formatStr = '%s'; +data = textscan(fID,formatStr,'Delimiter','\n'); +fclose(fID); fluxDataRaw = []; for i=1:length(data{1}) row = data{1}(i); row = strsplit(row{1},'\t'); row = row(1:end); + if i > 1 && numel(row) ~= size(fluxDataRaw,2) + error('loadFluxData:raggedRow', ... + 'Line %d of %s has %d cells, the header has %d.', ... + i, fluxDataFile, numel(row), size(fluxDataRaw,2)); + end fluxDataRaw = [fluxDataRaw; row]; end if size(fluxDataRaw,2) == 0 @@ -100,9 +109,15 @@ end %Extract observed byProduct names from file +if size(fluxDataRaw,2) < 4 + error('loadFluxData:noExchangeColumns', ... + '%s has no exchange flux columns after the columns condition, Ptot and grRate.', fluxDataFile); +end exchRxns = fluxDataRaw(1,4:end); -exchMets = strtrim(regexprep(exchRxns,'(.*)\(.*\)$','$1')); -exchRxns = regexprep(exchRxns,'.*\((.*)\)$','$1'); +exchMets = exchRxns; +for i=1:numel(exchRxns) + [exchMets{i}, exchRxns{i}] = splitExchHeader(exchRxns{i}); +end fluxData.conds = fluxDataRaw(2:end,1); fluxData.Ptot = str2double(fluxDataRaw(2:end,2)); @@ -112,3 +127,29 @@ fluxData.exchMets = exchMets; fluxData.exchRxnIDs = exchRxns; end + +function [met, rxn] = splitExchHeader(header) +% Split a column header ' ()' at the parenthesis that +% matches the final closing one, so a reaction id may contain parentheses. A +% header of another form gives the whole header as the reaction id. +met = strtrim(header); +rxn = header; +if ~endsWith(header,')') + return +end +depth = 0; +for k=numel(header):-1:1 + if header(k) == ')' + depth = depth + 1; + elseif header(k) == '(' + depth = depth - 1; + if depth == 0 + if k < numel(header)-1 + met = strtrim(header(1:k-1)); + rxn = header(k+1:end-1); + end + return + end + end +end +end diff --git a/src/kcat_tuning/evotune/loadEvotuneData.m b/src/kcat_tuning/evotune/loadEvotuneData.m index c2105b50..85c300ba 100644 --- a/src/kcat_tuning/evotune/loadEvotuneData.m +++ b/src/kcat_tuning/evotune/loadEvotuneData.m @@ -44,19 +44,26 @@ basePath = modelAdapter.params.path; -evoData.fluxData = loadFluxData(fullfile(basePath,'data','evotuneFluxData.tsv'), modelAdapter); +evoData.fluxData = []; +fluxFile = fullfile(basePath,'data','evotuneFluxData.tsv'); +if isfile(fluxFile) + evoData.fluxData = loadFluxData(fluxFile, modelAdapter); +end if ~isempty(evoData.fluxData) evoData.fluxData.biomass = modelAdapter.params.bioRxn; end -evoData.maxGrate = loadFluxData(fullfile(basePath,'data','evotuneMaxGrowth.tsv'), modelAdapter); +evoData.maxGrate = []; +maxGrateFile = fullfile(basePath,'data','evotuneMaxGrowth.tsv'); +if isfile(maxGrateFile) + evoData.maxGrate = loadFluxData(maxGrateFile, modelAdapter); +end if ~isempty(evoData.maxGrate) evoData.maxGrate.biomass = modelAdapter.params.bioRxn; end -% Unlike loadFluxData (used for fluxData/maxGrate above, which already -% degrades to [] on a missing file), readtable has no such guard, so this -% is wrapped explicitly to keep all three files uniformly optional. +% All three files are optional: a missing file gives [] (fluxData, maxGrate) +% or {} (zeroFlux). zeroExchFile = fullfile(basePath,'data','evotuneZeroExch.tsv'); if isfile(zeroExchFile) evoData.zeroFlux = table2cell(readtable(zeroExchFile, 'Delimiter', '\t', 'FileType','delimitedtext')); From 7d976fe3f2ee709ca0d099e276289581e7557a80 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 12/31] getEnzymeUsage leaves out enzymes without usage reaction; reportEnzymeUsage reads filtered usage and rejects a negative topAbsUsage --- src/enzyme_usage/getEnzymeUsage.m | 12 +++++++++++- src/enzyme_usage/reportEnzymeUsage.m | 15 +++++++++------ 2 files changed, 20 insertions(+), 7 deletions(-) diff --git a/src/enzyme_usage/getEnzymeUsage.m b/src/enzyme_usage/getEnzymeUsage.m index 32bd52be..ac09a085 100644 --- a/src/enzyme_usage/getEnzymeUsage.m +++ b/src/enzyme_usage/getEnzymeUsage.m @@ -39,7 +39,9 @@ % - capUsage : vector of enzyme capacity usages. % - absUsage : vector of absolute enzyme usages. % - UB : vector of enzyme exchange reaction upper bounds. -% - protID : string array of matching protein IDs. +% - protID : string array of matching protein IDs. An enzyme without a +% usage_prot_ reaction in the model is left out, with the warning +% getEnzymeUsage:missingUsageReaction. % - fluxes : vector of fluxes, copy of input fluxes. % % Examples @@ -60,6 +62,14 @@ end usageData.protID = ecModel.ec.enzymes; [~,rxnIdx] = ismember(strcat('usage_prot_',ecModel.ec.enzymes),ecModel.rxns); +noUsageRxn = rxnIdx == 0; +if any(noUsageRxn) + warning('getEnzymeUsage:missingUsageReaction', ... + 'The model has no usage reaction for the enzymes %s; they are left out.', ... + strjoin(usageData.protID(noUsageRxn), ', ')) + usageData.protID = usageData.protID(~noUsageRxn); + rxnIdx = rxnIdx(~noUsageRxn); +end usageData.UB = ecModel.ub(rxnIdx); usageData.absUsage = abs(fluxes(rxnIdx)); diff --git a/src/enzyme_usage/reportEnzymeUsage.m b/src/enzyme_usage/reportEnzymeUsage.m index 72cb88e5..ff313280 100644 --- a/src/enzyme_usage/reportEnzymeUsage.m +++ b/src/enzyme_usage/reportEnzymeUsage.m @@ -20,7 +20,7 @@ % how many top enzymes by absolute usage to list (default 10). Use 0 % or Inf to list all enzymes instead. Enzymes whose reactions carry no % flux are never listed here, as they do not contribute to the -% solution. +% solution. A negative value raises reportEnzymeUsage:invalidTopAbsUsage. % % Returns % ------- @@ -54,6 +54,9 @@ % values it used to reach the indexing below unclamped: usageData.protID(topUse(1:N)) % with N > numel(topUse) is an out-of-bounds index, so any model with fewer than 10 % enzymes (ecTestGEM has 5) crashed on the default call. +if ~(isnumeric(topAbsUsage) && isscalar(topAbsUsage) && topAbsUsage >= 0) + error('reportEnzymeUsage:invalidTopAbsUsage', 'topAbsUsage must be a non-negative number.') +end if topAbsUsage == 0 || isinf(topAbsUsage) || topAbsUsage > numel(usageData.protID) topAbsUsage = numel(usageData.protID); end @@ -84,8 +87,8 @@ if isscalar(find(carriedFlux)) highUsage.protID(end+1,1) = highEnzyme(i); highUsage.geneID(end+1,1) = geneIDs(i); - highUsage.absUsage(end+1,1) = usageData.absUsage(enzIdx(i)); - highUsage.capUsage(end+1,1) = usageData.capUsage(enzIdx(i)); + highUsage.absUsage(end+1,1) = usageData.absUsage(highUsageProt(i)); + highUsage.capUsage(end+1,1) = usageData.capUsage(highUsageProt(i)); highUsage.kcat(end+1,1) = kcat(carriedFlux); highUsage.source(end+1,1) = ecModel.ec.source(idx(carriedFlux)); highUsage.rxnID(end+1,1) = rxns(carriedFlux); @@ -95,8 +98,8 @@ % Add one entry for combined usage highUsage.protID(end+1,1) = highEnzyme(i); highUsage.geneID(end+1,1) = geneIDs(i); - highUsage.absUsage(end+1,1) = usageData.absUsage(enzIdx(i)); - highUsage.capUsage(end+1,1) = usageData.capUsage(enzIdx(i)); + highUsage.absUsage(end+1,1) = usageData.absUsage(highUsageProt(i)); + highUsage.capUsage(end+1,1) = usageData.capUsage(highUsageProt(i)); highUsage.kcat(end+1,1) = nan; highUsage.source{end+1,1} = '==='; highUsage.rxnID{end+1,1} = '==='; @@ -109,7 +112,7 @@ enzMet = strcat('prot_',highEnzyme{i}); [~, enzEcIdx] = ismember(enzMet,ecModel.mets); indAbsUse = full(transpose(-ecModel.S(enzEcIdx,rIdx)).*enzFlux); - indCapUse = (indAbsUse /sum(indAbsUse)) * usageData.capUsage(enzIdx(i)); + indCapUse = (indAbsUse /sum(indAbsUse)) * usageData.capUsage(highUsageProt(i)); rxnNumber = length(rIdx); highUsage.protID(end+1:end+rxnNumber,1) = highEnzyme(i); From 320118331a668b8d558f0533f4467d7b79a6cb16 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 13/31] truncateValues returns infinite values unchanged --- src/kcat_tuning/truncateValues.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/kcat_tuning/truncateValues.m b/src/kcat_tuning/truncateValues.m index f431baeb..f3fb3f47 100644 --- a/src/kcat_tuning/truncateValues.m +++ b/src/kcat_tuning/truncateValues.m @@ -2,7 +2,8 @@ % truncateValues Round values in selected columns to six significant digits. % % Rounds each value in the given columns of a cell array or table to six -% significant digits, based on the value's order of magnitude. +% significant digits, based on the value's order of magnitude. Infinite +% values are returned unchanged. % % Parameters % ---------- @@ -20,6 +21,9 @@ [m,n] = size(table); for i = 1:m for j = cols + if isinf(table{i,j}) + continue + end orderMagn = max([ceil(log10(abs(table{i,j}))),0]); table{i,j} = round(table{i,j},6-orderMagn); end From 6ee575a6df3de5a6bf96c2feb9b684f93174b5e7 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:06 +0100 Subject: [PATCH 14/31] Add tests for the contract changes of round 5 group 3 --- test/VerificationMatrix.tsv | 14 +- test/unit_tests/parityRound5G3Tests.m | 220 ++++++++++++++++++++++++++ 2 files changed, 231 insertions(+), 3 deletions(-) create mode 100644 test/unit_tests/parityRound5G3Tests.m diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index bf792183..5136f65f 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,6 +1,6 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time -change_model/applyComplexData.m tc0003, tc0004, tc0103 +change_model/applyComplexData.m tc0003, tc0004, tc0103, tc0501 change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238, tc0250 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. @@ -9,6 +9,8 @@ change_model/setProtPoolSize.m tc0005 change_model/getKcatAcrossIsoenzymes.m tc0011, tc0252 change_model/getSubsetEcModel.m tc0258, tc0259 change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 +enzyme_usage/getEnzymeUsage.m tc0507, tc0508 +enzyme_usage/reportEnzymeUsage.m tc0508, tc0509 gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m tc0256, tc0257 gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 @@ -33,11 +35,17 @@ gather_kcats/kcatsFromHomologs.m tc0095, tc0096 gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 change_model/upgradeEcModel.m tc0093, tc0094, tc0099 +io/loadFluxData.m tc0506 io/loadProtData.m tc0101, tc0102, tc0260 -change_model/addNewRxnsToEC.m tc0211, tc0232, tc0237 +change_model/addNewRxnsToEC.m tc0211, tc0232, tc0237, tc0500 change_model/mapRxnsToConv.m tc0213 flux_analysis/ecFVA.m tc0214 flux_analysis/ecFSEOF.m tc0100, tc0236 +limit_proteins/calculateFfactor.m tc0502 +limit_proteins/constrainEnzConcs.m tc0503 +limit_proteins/constrainFluxData.m tc0505 +limit_proteins/fillEnzConcs.m tc0504 +limit_proteins/relaxProteomicsGreedy.m tc0511 utilities/upgradeGECKOproject.m tc0094, tc0099 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases @@ -57,7 +65,7 @@ kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed fo kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253, tc0254 kcat_tuning/sigmaFitter.m tc0020, tc0255 kcat_tuning/topUsedEnzymes.m Not explicitly tested - The code has existed for a long time Not sure this works anymore - written for the old model structure? -kcat_tuning/truncateValues.m Not explicitly tested - The code has existed for a long time +kcat_tuning/truncateValues.m tc0510 limit_proteins/constrainProtConcs.m tc0013 limit_proteins/findLimitingUBs.m Not explicitly tested - The code has existed for a long time Probably not relevant anymore, remove? limit_proteins/fitGAM.m diff --git a/test/unit_tests/parityRound5G3Tests.m b/test/unit_tests/parityRound5G3Tests.m new file mode 100644 index 00000000..a9097005 --- /dev/null +++ b/test/unit_tests/parityRound5G3Tests.m @@ -0,0 +1,220 @@ +%run this test case with the command +%results = runtests('parityRound5G3Tests.m') + +function tests = parityRound5G3Tests + tests = functiontests(localfunctions); +end + +function testAddNewRxnsToECWithoutNewEnzymes_tc0500(testCase) + % The genes of the new reaction already have enzymes: the reaction is + % added and no enzyme is. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + nEnzymes = numel(ecModel.ec.enzymes); + nEcRxns = numel(ecModel.ec.rxns); + newRxns.rxns = {'RNEW'}; + newRxns.rxnNames = {'RNEW'}; + newRxns.equations = {'m1[c] => e2[e]'}; + newRxns.grRules = {'G1'}; + [ecModel, rxnsAdded, enzAdded] = addNewRxnsToEC(ecModel, newRxns, [], adapter); + verifyEqual(testCase, rxnsAdded, {'RNEW'}) + verifyEmpty(testCase, enzAdded) + verifyEqual(testCase, numel(ecModel.ec.enzymes), nEnzymes) + verifyEqual(testCase, numel(ecModel.ec.rxns), nEcRxns + 1) + [~, rxnIdx] = ismember('RNEW', ecModel.ec.rxns); + [~, enzIdx] = ismember('P1', ecModel.ec.enzymes); + verifyEqual(testCase, full(ecModel.ec.rxnEnzMat(rxnIdx, :)), double((1:nEnzymes) == enzIdx)) +end + +function testApplyComplexDataPartialZeroStoichiometry_tc0501(testCase) + % A subunit with stoichiometry 0 counts as a subunit with stoichiometry 1. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + complexInfo = jsondecode(['[{"complexID":"R2Compl","name":"R2Compl","specie":"x",' ... + '"geneName":["G1","G2"],"protID":["P1","P2"],"stochiometry":[2,0]}]']); + [ecModel, foundComplex] = applyComplexData(ecModel, 'complexInfo', complexInfo, ... + 'modelAdapter', adapter, 'verbose', false); + [~, p12] = ismember({'P1','P2'}, ecModel.ec.enzymes); + verifyEqual(testCase, full(ecModel.ec.rxnEnzMat(1, p12)), [2 1]) + verifyEqual(testCase, height(foundComplex), 2) +end + +function testCalculateFfactorEmptyEnzymeList_tc0502(testCase) + % An empty cell array means no enzymes, [] means every enzyme of the model. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + protData.uniprotIDs = {'P1';'P2';'X1'}; + protData.abundances = [1;1;2]; + verifyEqual(testCase, calculateFfactor(ecModel, 'protData', protData, ... + 'enzymes', {}, 'modelAdapter', adapter), 0) + verifyEqual(testCase, calculateFfactor(ecModel, 'protData', protData, ... + 'modelAdapter', adapter), 0.5) + verifyEqual(testCase, calculateFfactor(ecModel, 'protData', protData, ... + 'enzymes', {'P1'}, 'modelAdapter', adapter), 0.25) +end + +function testConstrainEnzConcsNegativeConcentrationRaises_tc0503(testCase) + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel.ec.concs(:) = NaN; + ecModel.ec.concs(1) = -0.5; + verifyError(testCase, @() constrainEnzConcs(ecModel), 'constrainEnzConcs:negativeConcentration') + % removing the constraints does not read the concentrations + cleared = constrainEnzConcs(ecModel, 'removeConstraints', true); + verifyEqual(testCase, cleared.ub(startsWith(cleared.rxns, 'usage_prot_')), ... + 1000 * ones(numel(ecModel.ec.enzymes), 1)) +end + +function testFillEnzConcsNegativeAbundanceRaises_tc0504(testCase) + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + protData.uniprotIDs = {'P1';'X1'}; + protData.abundances = [-0.5;-1]; + verifyError(testCase, @() fillEnzConcs(ecModel, protData), 'fillEnzConcs:negativeAbundance') + % a negative abundance of a protein that is not an enzyme of the model is not read + protData.uniprotIDs = {'P1';'X1'}; + protData.abundances = [0.5;-1]; + ecModel = fillEnzConcs(ecModel, protData); + verifyEqual(testCase, ecModel.ec.concs(1), 0.5) +end + +function testConstrainFluxDataInvalidArguments_tc0505(testCase) + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + fluxData.conds = {'c1'}; + fluxData.grRate = 0.1; + fluxData.exchFluxes = 1; + fluxData.exchRxnIDs = {'S1'}; + fluxData.exchMets = {'m'}; + constrained = constrainFluxData(ecModel, 'fluxData', fluxData, 'modelAdapter', adapter); + verifyEqual(testCase, constrained.ub(strcmp(constrained.rxns, 'S1')), 1) + verifyError(testCase, @() constrainFluxData(ecModel, 'fluxData', fluxData, ... + 'maxMinGrowth', 'median', 'modelAdapter', adapter), 'constrainFluxData:invalidMaxMinGrowth') + fluxData.exchRxnIDs = {'NOT_A_RXN'}; + verifyError(testCase, @() constrainFluxData(ecModel, 'fluxData', fluxData, ... + 'modelAdapter', adapter), 'constrainFluxData:missingReaction') + fluxData.exchRxnIDs = {'S1'}; + verifyError(testCase, @() constrainFluxData(ecModel, 'fluxData', fluxData, ... + 'condition', 'nope', 'modelAdapter', adapter), 'constrainFluxData:unknownCondition') +end + +function testLoadFluxDataMalformedFiles_tc0506(testCase) + adapter = getAdapter(); + folder = tempname; + mkdir(folder) + cleanup = onCleanup(@() rmdir(folder, 's')); + tab = char(9); + + verifyError(testCase, @() loadFluxData(fullfile(folder, 'missing.tsv'), adapter), ... + 'loadFluxData:fileNotFound') + + fixedOnly = writeLines(fullfile(folder, 'fixed.tsv'), ... + {['cond' tab 'Ptot' tab 'grRate'], ['c1' tab '0.5' tab '0.4']}); + verifyError(testCase, @() loadFluxData(fixedOnly, adapter), 'loadFluxData:noExchangeColumns') + + ragged = writeLines(fullfile(folder, 'ragged.tsv'), ... + {['cond' tab 'Ptot' tab 'grRate' tab 'Glc (EX_a)' tab 'O2 (EX_b)'], ... + ['c1' tab '0.5' tab '0.4' tab '-1' tab '2'], ... + ['c2' tab '0.5' tab '0.4' tab '-1']}); + verifyError(testCase, @() loadFluxData(ragged, adapter), 'loadFluxData:raggedRow') + + nested = writeLines(fullfile(folder, 'nested.tsv'), ... + {['cond' tab 'Ptot' tab 'grRate' tab 'X (EX_a(e))' tab 'O2 (EX_b)'], ... + ['c1' tab '0.5' tab '0.4' tab '-1' tab '2']}); + fluxData = loadFluxData(nested, adapter); + verifyEqual(testCase, fluxData.exchMets, {'X', 'O2'}) + verifyEqual(testCase, fluxData.exchRxnIDs, {'EX_a(e)', 'EX_b'}) +end + +function testGetEnzymeUsageLeavesOutMissingUsageReaction_tc0507(testCase) + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = removeReactions(ecModel, 'usage_prot_P3'); + fluxes = zeros(numel(ecModel.rxns), 1); + verifyWarning(testCase, @() getEnzymeUsage(ecModel, fluxes), 'getEnzymeUsage:missingUsageReaction') + warning('off', 'getEnzymeUsage:missingUsageReaction') + restore = onCleanup(@() warning('on', 'getEnzymeUsage:missingUsageReaction')); + usageData = getEnzymeUsage(ecModel, fluxes); + verifyEqual(testCase, usageData.protID, setdiff(ecModel.ec.enzymes, {'P3'}, 'stable')) + verifyEqual(testCase, numel(usageData.absUsage), numel(usageData.protID)) +end + +function testReportEnzymeUsageWithoutZeroUsage_tc0508(testCase) + % getEnzymeUsage with zero=false filters the enzymes; the report reads + % the usage of the filtered enzymes. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + ecModel.ub(strcmp(ecModel.rxns, 'usage_prot_P5')) = 10; + fluxes = zeros(numel(ecModel.rxns), 1); + fluxes(strcmp(ecModel.rxns, 'usage_prot_P5')) = 10; + fluxes(strcmp(ecModel.rxns, 'R5')) = 10; + usageData = getEnzymeUsage(ecModel, fluxes, 'zero', false); + verifyEqual(testCase, usageData.protID, {'P5'}) + report = reportEnzymeUsage(ecModel, usageData); + verifyEqual(testCase, report.highCapUsage.protID, {'P5'}) + verifyEqual(testCase, report.highCapUsage.absUsage, 10) + verifyEqual(testCase, report.highCapUsage.capUsage, 1) + verifyEqual(testCase, report.topAbsUsage.protID, {'P5'}) +end + +function testReportEnzymeUsageNegativeTopAbsUsageRaises_tc0509(testCase) + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 100; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, [], [], [], adapter); + usageData = getEnzymeUsage(ecModel, zeros(numel(ecModel.rxns), 1)); + verifyError(testCase, @() reportEnzymeUsage(ecModel, usageData, 'topAbsUsage', -2), ... + 'reportEnzymeUsage:invalidTopAbsUsage') +end + +function testTruncateValuesInfinite_tc0510(testCase) + values = {Inf, -Inf, 1234565, 0.5}; + truncated = truncateValues(values, 1:4); + verifyEqual(testCase, truncated{1}, Inf) + verifyEqual(testCase, truncated{2}, -Inf) + verifyEqual(testCase, truncated{3}, 1234570) + verifyEqual(testCase, truncated{4}, 0.5) +end + +function testRelaxProteomicsGreedyEmptyEnzymeSet_tc0511(testCase) + % An empty cell array makes no enzyme eligible; [] makes every + % constrained enzyme eligible. + adapter = getAdapter(); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 10; + ecModel.ec.source(:) = {'manual'}; + ecModel = applyKcatConstraints(ecModel); + ecModel = setProtPoolSize(ecModel, 0.5, 0.5, 0.5, adapter); + ecModel.ec.concs(:) = NaN; + ecModel.ec.concs(strcmp(ecModel.ec.enzymes, 'P5')) = 10; + ecModel = constrainEnzConcs(ecModel); + + warning('off', 'all') + restore = onCleanup(@() warning('on', 'all')); + none = relaxProteomicsGreedy(ecModel, 'minimalGrowth', 50, 'enzymeSet', {}); + verifyFalse(testCase, none.converged) + verifyEmpty(testCase, fieldnames(none.relaxed)) + all = relaxProteomicsGreedy(ecModel, 'minimalGrowth', 50, 'enzymeSet', []); + verifyTrue(testCase, all.converged) + verifyEqual(testCase, fieldnames(all.relaxed), {'P5'}) +end + +function adapter = getAdapter() +geckoPath = findGECKOroot; +adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function file = writeLines(file, lines) +fid = fopen(file, 'w'); +fprintf(fid, '%s\n', lines{:}); +fclose(fid); +end From 78e17572ed1e90d04cedd50c90f176c3739c6256 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:23 +0100 Subject: [PATCH 15/31] getECfromDatabase treats an empty logical ecRxns as no selection --- src/get_enzyme_data/getECfromDatabase.m | 10 ++++++++-- test/VerificationMatrix.tsv | 2 +- test/unit_tests/ecFromDatabaseConflictTests.m | 15 +++++++++++++++ 3 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/get_enzyme_data/getECfromDatabase.m b/src/get_enzyme_data/getECfromDatabase.m index 16e88f3f..102fe838 100644 --- a/src/get_enzyme_data/getECfromDatabase.m +++ b/src/get_enzyme_data/getECfromDatabase.m @@ -16,8 +16,10 @@ % of length model.ec.rxns that specifies which model.ec.eccodes entries % should be queried. Existing values in model.ec.eccodes will be wiped. % Entries that are indicated by false will be kept and not modified by -% this function, and are not searched for multiplicities (by default all -% model.ec.eccodes entries are populated by this function). +% this function, and are not searched for multiplicities. An empty +% logical array selects no entry, and the model is returned unchanged. +% When ecRxns is omitted or [], all model.ec.eccodes entries are +% populated by this function. % action : char % response action if multiple proteins with different EC numbers are % found for a given gene in a metabolic reaction (default 'display'): @@ -62,6 +64,10 @@ end params = modelAdapter.getParameters(); +if islogical(ecRxns) && isempty(ecRxns) + return +end + rxnEnzMat = model.ec.rxnEnzMat; genes = modelAdapter.getUniprotCompatibleGenes(model.ec.genes); diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 2ea31f48..3244566f 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -45,7 +45,7 @@ gather_kcats/ecStatsFromBrenda.m tc0400, tc0401 - not verified against a real BR gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243, tc0244, tc0245, tc0246 -get_enzyme_data/getECfromDatabase.m tc0007, tc0404, tc0405 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM +get_enzyme_data/getECfromDatabase.m tc0007, tc0404, tc0405, tc0406 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested; also tc0246 download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 diff --git a/test/unit_tests/ecFromDatabaseConflictTests.m b/test/unit_tests/ecFromDatabaseConflictTests.m index a8cc6027..a511106f 100644 --- a/test/unit_tests/ecFromDatabaseConflictTests.m +++ b/test/unit_tests/ecFromDatabaseConflictTests.m @@ -26,6 +26,21 @@ function testGetECfromDatabaseIgnoresConflictOutsideSelection_tc0405(testCase) verifyEqual(testCase, ecModel2.ec.eccodes(~selected), repmat({''}, nnz(~selected), 1)) end +function testGetECfromDatabaseEmptySelectionChangesNothing_tc0406(testCase) + % An empty logical selection searches nothing and returns the model as it + % is; [] selects every ec reaction. + [ecModel, adapter] = conflictModel(testCase); + ecModel.ec.eccodes = repmat({'9.9.9.9'}, numel(ecModel.ec.rxns), 1); + lastwarn('') + same = getECfromDatabase(ecModel, false(0,1), 'display', adapter); + verifyEqual(testCase, same, ecModel) + verifyEmpty(testCase, lastwarn()) + warning('off', 'getECfromDatabase:conflicts') + cleanup = onCleanup(@() warning('on', 'getECfromDatabase:conflicts')); %#ok + all = getECfromDatabase(ecModel, [], 'ignore', adapter); + verifyEqual(testCase, all.ec.eccodes{1}, '1.1.1.1') +end + function [ecModel, adapter] = conflictModel(testCase) % ecTestGEM with a second UniProt entry for gene G4 and another EC number. geckoPath = findGECKOroot; From 8ab8bd7e6d78a51ac07b8c61a713876967e89649 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:11 +0100 Subject: [PATCH 16/31] collectKcats counts each source file once, and mergeKcats matches OKP-tagged sources by their underlying source --- src/gather_kcats/collectKcats.m | 8 +++-- src/gather_kcats/mergeKcats.m | 8 +++-- test/VerificationMatrix.tsv | 4 +-- test/unit_tests/kcatContractTests.m | 55 +++++++++++++++++++++++++++++ 4 files changed, 69 insertions(+), 6 deletions(-) create mode 100644 test/unit_tests/kcatContractTests.m diff --git a/src/gather_kcats/collectKcats.m b/src/gather_kcats/collectKcats.m index 7ae960f3..2bf88084 100644 --- a/src/gather_kcats/collectKcats.m +++ b/src/gather_kcats/collectKcats.m @@ -198,6 +198,7 @@ continue end lists = {}; + usedFiles = {}; priority = {}; for t = 1:numel(st.files) token = st.files{t}; @@ -209,8 +210,9 @@ if isempty(list) continue end - if ~any(cellfun(@(l) isequal(l, list), lists)) + if ~ismember(name, usedFiles) lists{end+1} = list; %#ok + usedFiles{end+1} = name; %#ok end if startsWith(token, 'database_') priority{end+1} = token; %#ok @@ -284,10 +286,12 @@ end function s = normalizeSource(label) -% As mergeKcats folds source labels: lower case, other characters as '_'. +% As mergeKcats folds source labels: lower case, other characters as '_', +% without a leading okp_ (the tag of OpenKineticsPredictor rows). s = lower(strtrim(char(label))); s = regexprep(s, '[^0-9a-z]+', '_'); s = regexprep(s, '^_+|_+$', ''); +s = regexprep(s, '^okp_', ''); end function name = sourceFile(token) diff --git a/src/gather_kcats/mergeKcats.m b/src/gather_kcats/mergeKcats.m index 483c51f3..4930b305 100644 --- a/src/gather_kcats/mergeKcats.m +++ b/src/gather_kcats/mergeKcats.m @@ -27,8 +27,11 @@ % wildcardLimit / bottomOriginLimit. % % Any other token is matched against each row's source after folding -% to lowercase snake_case (e.g. 'dlkcat', 'catapro'). Sources not -% listed are dropped. +% to lowercase snake_case (e.g. 'dlkcat', 'catapro'). A leading 'okp_' +% of the row's folded source (the OKP- tag of OpenKineticsPredictor +% rows) is dropped for matching, so 'OKP-BRENDA' routes to the database +% tiers and 'OKP-CataPro' matches 'catapro'; the rows keep their +% source label. Sources not listed are dropped. % % Name-Value Arguments % -------------------- @@ -170,6 +173,7 @@ end token = cellfun(@normalizeSource, rawSource, 'UniformOutput', false); +token = regexprep(token, '^okp_', ''); isDb = ismember(token, dbSources); hasMeta = ~isnan(wcl) & ~isnan(origin); diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index bf792183..c9ad2ac8 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -22,7 +22,7 @@ gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 gather_kcats/kcatsFromBrenda.m tc0074, tc0078, tc0226 gather_kcats/kcatsFromDLKcat.m tc0074, tc0078 - reads an existing DLKcat output; running DLKcat needs Docker gather_kcats/kcatsFromOKP.m tc0075 - reads a stored result; submitting needs the OKP service -gather_kcats/collectKcats.m tc0079, tc0080, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090, tc0098 +gather_kcats/collectKcats.m tc0079, tc0080, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090, tc0098, tc0300 gather_kcats/assignKcats.m tc0079, tc0083 gather_kcats/updateKcatSources.m tc0082, tc0087 gather_kcats/kcatAssignmentConfig.m tc0079, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090 @@ -42,7 +42,7 @@ utilities/upgradeGECKOproject.m tc0094, tc0099 gather_kcats/readOpenKineticsPredictorOutput.m tc0030, tc0069 gather_kcats/reviewKcatAssignment.m Has only been tested by running it manually, not tested by test cases gather_kcats/ecStatsFromBrenda.m Not verified against real BRENDA data in this environment -- see its own TODO comment -gather_kcats/mergeKcats.m tc0240 +gather_kcats/mergeKcats.m tc0240, tc0301 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243, tc0244, tc0245, tc0246 get_enzyme_data/getECfromDatabase.m tc0007 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM diff --git a/test/unit_tests/kcatContractTests.m b/test/unit_tests/kcatContractTests.m new file mode 100644 index 00000000..c90fa722 --- /dev/null +++ b/test/unit_tests/kcatContractTests.m @@ -0,0 +1,55 @@ +%run this test case with the command +%results = runtests('kcatContractTests.m') + +function tests = kcatContractTests + tests = functiontests(localfunctions); +end + +function testCollectKcatsCountsEachSourceFileOnce_tc0300(testCase) + % database_top and database_bottom both read brenda.tsv, which has an + % empty log10_sd (NaN) column: the file is one list, so R3 has the two + % values of the file as candidates, not four. + [model, adapter, project] = kcatProject(); + adapter.params.kcat_assignment = struct('isozymes', false, 'standard', false, 'stage', ... + struct('action', 'select', 'sources', {{'database_top', 'database_bottom'}})); + [baseRxns, enz] = getEcRxnKeys(model); + i = find(strcmp(model.ec.rxns, 'R3')); + writeLines(fullfile(project, 'data', 'kcats', 'brenda.tsv'), {'# source: brenda', ... + sprintf('rxn\tenzymes\tsubstrate\tkcat\tlog10_sd\tsource\tdetail'), ... + sprintf('%s\t%s\tm1\t5\t\tbrenda\twc=0; origin=1', baseRxns{i}, enz{i}), ... + sprintf('%s\t%s\tm1\t20\t\tbrenda\twc=0; origin=1', baseRxns{i}, enz{i})}); + decisions = collectKcats(model, 'modelAdapter', adapter, 'writeReport', false); + row = strcmp(decisions.ecRxn, 'R3'); + verifyEqual(testCase, decisions.kcat(row), 20) + verifyEqual(testCase, decisions.candidates(row), 2) + verifyEqual(testCase, decisions.spreadLog10(row), log10(4), 'AbsTol', 1e-12) +end + +function testMergeKcatsMatchesOkpTaggedSources_tc0301(testCase) + % The OKP- tag is dropped for matching: OKP-BRENDA is an exact database + % hit and OKP-CataPro matches the priority entry catapro; the rows keep + % their labels and win over DLKcat. + list = struct('source', 'x', 'rxns', {{'A'; 'A'; 'B'; 'B'}}, 'kcats', [3; 9; 4; 9], ... + 'kcatSource', {{'OKP-BRENDA'; 'DLKcat'; 'OKP-CataPro'; 'DLKcat'}}); + merged = mergeKcats({list}, {'database_exact', 'catapro', 'dlkcat'}); + verifyEqual(testCase, merged.rxns, {'A'; 'B'}) + verifyEqual(testCase, merged.kcats, [3; 4]) + verifyEqual(testCase, merged.kcatSource, {'OKP-BRENDA'; 'OKP-CataPro'}) +end + +function [model, adapter, project] = kcatProject() + adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); + project = tempname; + copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); + delete(fullfile(project, 'data', 'customKcats.tsv')); + mkdir(fullfile(project, 'data', 'kcats')); + adapter.params.path = project; + model = makeEcModel(getGeckoTestModel(), false, adapter); + model = getECfromGEM(model); +end + +function writeLines(file, lines) + fid = fopen(file, 'w'); + fprintf(fid, '%s\n', lines{:}); + fclose(fid); +end From df8e34c4edc07f9c4e4b1d734ee214bab64a6702 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:11 +0100 Subject: [PATCH 17/31] loadKcatSource skips a row of empty quoted fields and raises for a field with an opening quote and no closing one --- src/gather_kcats/loadKcatSource.m | 18 +++++++++++++++--- test/VerificationMatrix.tsv | 2 +- test/unit_tests/kcatContractTests.m | 19 +++++++++++++++++++ 3 files changed, 35 insertions(+), 4 deletions(-) diff --git a/src/gather_kcats/loadKcatSource.m b/src/gather_kcats/loadKcatSource.m index 6169d23f..105462f3 100644 --- a/src/gather_kcats/loadKcatSource.m +++ b/src/gather_kcats/loadKcatSource.m @@ -21,8 +21,10 @@ % % The file is read as UTF-8 (a leading byte-order mark is ignored). The % '#' lines before the column header are the header; '#' lines and blank -% lines further down are skipped. Fields wrapped in double quotes, as a -% spreadsheet may save them, are unwrapped. kcat and log10_sd must be plain +% lines further down are skipped, as are rows without a filled field (a row +% of empty quoted fields included). Fields wrapped in double quotes, as a +% spreadsheet may save them, are unwrapped; a field that opens a double quote +% and does not close it is an error. kcat and log10_sd must be plain % numbers (12, 0.5, 1e-05); other text, such as '5,5' with a decimal comma, % is reported as invalid. A row with more filled fields than there are % columns is an error. @@ -263,7 +265,11 @@ nCol = numel(columns); trimmed = strtrim(lines); -skip = cellfun(@isempty, trimmed) | startsWith(trimmed, '#') | startsWith(trimmed, '"#'); +skip = cellfun(@isempty, trimmed) | startsWith(trimmed, '#'); +commentLike = find(startsWith(trimmed, '"#')); +for q = commentLike(:)' + skip(q) = startsWith(unquoteField(trimmed{q}), '#'); +end skip(1:k) = true; kept = find(~skip); L = string(lines(kept)); @@ -289,6 +295,12 @@ end F = strtrim(F); quoted = find(startsWith(F, '"')); +unclosed = quoted(~(strlength(F(quoted)) >= 2 & endsWith(F(quoted), '"'))); +if ~isempty(unclosed) + badRow = min(mod(unclosed - 1, size(F, 1)) + 1); + error('loadKcatSource:quote', '%s, line %d has a field with an opening double quote and no closing one.', ... + filename, kept(badRow)) +end for q = quoted(:)' F(q) = unquoteField(char(F(q))); end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index c9ad2ac8..1aec145e 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -18,7 +18,7 @@ gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameteri gather_kcats/writeDLKcatInput.m partly by tc0011 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 -gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078 +gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078, tc0302 gather_kcats/kcatsFromBrenda.m tc0074, tc0078, tc0226 gather_kcats/kcatsFromDLKcat.m tc0074, tc0078 - reads an existing DLKcat output; running DLKcat needs Docker gather_kcats/kcatsFromOKP.m tc0075 - reads a stored result; submitting needs the OKP service diff --git a/test/unit_tests/kcatContractTests.m b/test/unit_tests/kcatContractTests.m index c90fa722..641353d4 100644 --- a/test/unit_tests/kcatContractTests.m +++ b/test/unit_tests/kcatContractTests.m @@ -37,6 +37,25 @@ function testMergeKcatsMatchesOkpTaggedSources_tc0301(testCase) verifyEqual(testCase, merged.kcatSource, {'OKP-BRENDA'; 'OKP-CataPro'}) end +function testLoadKcatSourceQuotedRows_tc0302(testCase) + % A row of empty quoted fields is blank and a quoted comment is a + % comment; a field that opens a quote and does not close it is an error. + [model, ~, project] = kcatProject(); + file = fullfile(project, 'odd.tsv'); + writeLines(file, {'# source: odd', sprintf('rxn\tenzymes\tkcat\tsource'), ... + sprintf('R3\tP4\t25\todd'), '"# a quoted comment"', sprintf('""\t""\t""\t""'), ... + sprintf('R5\tP5\t30\todd')}); + [list, report] = loadKcatSource(model, file); + verifyEqual(testCase, numel(report.records.kcat), 2) + verifyEmpty(testCase, report.invalid) + writeLines(file, {'# source: odd', sprintf('rxn\tenzymes\tkcat\tsource'), ... + sprintf('R3\tP4\t25\todd'), '"# a comment that was quoted'}); + verifyError(testCase, @() loadKcatSource(model, file), 'loadKcatSource:quote') + writeLines(file, {'# source: odd', sprintf('rxn\tenzymes\tkcat\tsource'), ... + sprintf('"R3\tP4\t25\todd')}); + verifyError(testCase, @() loadKcatSource(model, file), 'loadKcatSource:quote') +end + function [model, adapter, project] = kcatProject() adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); project = tempname; From 406c18e68087e12d23f1334888ca2f7766fb21c6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:11 +0100 Subject: [PATCH 18/31] readDLKcatOutput drops rows without a finite kcat before checking reactions and leaves out unknown substrates; writeDLKcatInput writes an empty file for an empty selection --- src/gather_kcats/readDLKcatOutput.m | 55 ++++++++++++++------- src/gather_kcats/writeDLKcatInput.m | 25 +++++++--- test/VerificationMatrix.tsv | 4 +- test/unit_tests/geckoCoreFunctionTests.m | 2 +- test/unit_tests/kcatContractTests.m | 62 ++++++++++++++++++++++++ 5 files changed, 122 insertions(+), 26 deletions(-) diff --git a/src/gather_kcats/readDLKcatOutput.m b/src/gather_kcats/readDLKcatOutput.m index f1cadc9f..d9aff5d3 100644 --- a/src/gather_kcats/readDLKcatOutput.m +++ b/src/gather_kcats/readDLKcatOutput.m @@ -16,6 +16,15 @@ % the obj.params.path folder specified in the modelAdapter). % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). +% strict : logical +% raise an error for a substrate that is not in model.metNames +% (readDLKcatOutput:unknownSubstrate). By default the rows with such a +% substrate are left out with a warning (default false). +% +% Rows without a numeric kcat ('NA', 'None', empty, text, 'inf') are left out +% before the reactions are checked. A file without any such row raises +% readDLKcatOutput:noKcats. A reaction that is not in model.ec.rxns raises +% readDLKcatOutput:unknownReaction. % % Returns % ------- @@ -45,9 +54,11 @@ p = parseGECKOargs(varargin, { ... 'outFile', []; ... - 'modelAdapter', []}); + 'modelAdapter', []; ... + 'strict', false}); outFile = p.outFile; modelAdapter = p.modelAdapter; +strict = p.strict; if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -65,36 +76,46 @@ DLKcatOutput = textscan(fID,'%s %s %s %s %s %s','Delimiter','\t','HeaderLines',1); fclose(fID); -% Check that DLKcat output file and model match (not fool proof, but good enough) [rxns, genes, subs, kcats] = deal(DLKcatOutput{[1,2,3,6]}); +% A row has a kcat when its kcat field is a plain finite number. +isNumber = ~cellfun(@isempty, regexp(strtrim(kcats), '^[-+]?(\d+\.?\d*|\.\d+)([eE][-+]?\d+)?$', 'once')); +kcatValues = nan(numel(kcats), 1); +kcatValues(isNumber) = str2double(kcats(isNumber)); +usable = isfinite(kcatValues); + % Check if it contains any kcat values -if all(cellfun(@isempty,kcats)) || all(strcmpi(kcats,'NA')) - error('DLKcat file does not contain any kcat values, please run runDLKcat() first.') +if ~any(usable) + error('readDLKcatOutput:noKcats', 'DLKcat file does not contain any kcat values, please run runDLKcat() first.') end -% Check that all substrates are in the model (case-insensitive: an SBML loader can +% Check that the substrates are in the model (case-insensitive: an SBML loader can % normalise metabolite-name capitalisation differently than whatever produced the -% DLKcat input, and that must not read as a genuinely missing substrate). +% DLKcat input, and that must not read as a genuinely missing substrate). Rows with +% another substrate are left out with a warning, or raise an error with 'strict'. matchMets = ismember(lower(subs),lower(model.metNames)); if ~all(matchMets) errorText = 'DLKcat was likely run with an input file that was generated from another ecModel, as the following substrates from DLKcat output cannot be found in model.metNames:'; - error('%s', ravenList(errorText,subs(~matchMets),false)) + if strict + error('readDLKcatOutput:unknownSubstrate', '%s', ravenList(errorText,subs(~matchMets),false)) + end + warning('readDLKcatOutput:unknownSubstrate', 'Rows with a substrate that is not in model.metNames are left out: %s', ... + strjoin(unique(subs(~matchMets))', ', ')) + usable = usable & matchMets; end -% Check that all reactions are in model.ec.rxns +% Check that the reactions of the rows with a kcat are in model.ec.rxns matchRxns = ismember(rxns,model.ec.rxns); -if ~all(matchRxns) - errorText = 'DLKcat was likely run with an input file that was generated from another ecModel, as the following reactions from DLKcat output cannot be found in model.metNames:'; - error('%s', ravenList(errorText,rxns(~matchRxns),false)) +if ~all(matchRxns(usable)) + errorText = 'DLKcat was likely run with an input file that was generated from another ecModel, as the following reactions from DLKcat output cannot be found in model.ec.rxns:'; + error('readDLKcatOutput:unknownReaction', '%s', ravenList(errorText,rxns(usable & ~matchRxns),false)) end -% Filter out entries with no numeric value -noOutput = cellfun(@isempty,regexpi(kcats,'[0-9]')); -kcats = str2double(kcats(~noOutput)); -rxns(noOutput) = []; -genes(noOutput) = []; -subs(noOutput) = []; +% Keep the entries with a numeric value +kcats = kcatValues(usable); +rxns = rxns(usable); +genes = genes(usable); +subs = subs(usable); % Make kcatList structure kcatList.source = 'DLKcat'; diff --git a/src/gather_kcats/writeDLKcatInput.m b/src/gather_kcats/writeDLKcatInput.m index ee2b9314..1760a576 100644 --- a/src/gather_kcats/writeDLKcatInput.m +++ b/src/gather_kcats/writeDLKcatInput.m @@ -14,7 +14,9 @@ % ecRxns : logical % for which reactions (from model.ec.rxns) DLKcat should predict kcat % values, provided as logical vector with same length as model.ec.rxns -% (default: all reactions). +% (default: all reactions). A vector without true entries selects no +% reaction and gives an empty file. The rows follow the order of +% model.ec.rxns. % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % onlyWithSmiles : logical @@ -56,7 +58,7 @@ if isempty(ecRxnsSel) ecRxnsSel = true(numel(model.ec.rxns),1); -elseif ~logical(ecRxnsSel) +elseif ~islogical(ecRxnsSel) error('ecRxns should be provided as logical vector') elseif numel(ecRxnsSel)~=numel(model.ec.rxns) error('Length of ecRxns is not the same as model.ec.rxns') @@ -91,6 +93,12 @@ 'modelAdapter', modelAdapter); % Prepare output +if isempty(pairs.ecRxn) + out = cell(6,0); + writeRows(filename, out) + writtenTable = out; + return +end out(1,:) = model.ec.rxns(pairs.ecRxn); out(2,:) = model.ec.genes(pairs.protein); out(3,:) = model.metNames(pairs.met); @@ -109,11 +117,16 @@ out(6,:) = cell(numel(out(1,:)),1); out(6,:) = {'NA'}; -% Write file +writeRows(filename, out) +writtenTable = out; +end + +function writeRows(filename, out) +% Write the rows (one per column of out); an empty table gives an empty file. fID = fopen(filename,'w'); -fprintf(fID,'%s\t%s\t%s\t%s\t%s\t%s\n',out{:}); +if ~isempty(out) + fprintf(fID,'%s\t%s\t%s\t%s\t%s\t%s\n',out{:}); +end fclose(fID); fprintf('Model-specific DLKcat input stored at %s\n',filename); - -writtenTable = out; end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 1aec145e..22111205 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -12,10 +12,10 @@ change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 gather_kcats/getStandardKcat.m tc0256, tc0257 gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 -gather_kcats/readDLKcatOutput.m tcm0001, tcm0002 - The test doesn't really check the function, only that some output is read +gather_kcats/readDLKcatOutput.m tcm0001, tcm0002, tc0303, tc0304 - The test doesn't really check the function, only that some output is read gather_kcats/runDLKcat.m tcm0001, tcm0002 - The test doesn't really check the function, only that we get something gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameterizations here; also tc0104 -gather_kcats/writeDLKcatInput.m partly by tc0011 +gather_kcats/writeDLKcatInput.m partly by tc0011, tc0305 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 gather_kcats/loadKcatSource.m tc0071, tc0072, tc0073, tc0076, tc0077, tc0078, tc0302 diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 308b830f..d4962dc0 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -1422,7 +1422,7 @@ function testReadDLKcatOutputUnrecognizedSubstrateErrorsCleanly_tc0037(testCase) fclose(fid); try - readDLKcatOutput(ecModel, 'outFile', outFile, 'modelAdapter', adapter); + readDLKcatOutput(ecModel, 'outFile', outFile, 'modelAdapter', adapter, 'strict', true); errored = false; msg = ''; catch ME diff --git a/test/unit_tests/kcatContractTests.m b/test/unit_tests/kcatContractTests.m index 641353d4..32f19998 100644 --- a/test/unit_tests/kcatContractTests.m +++ b/test/unit_tests/kcatContractTests.m @@ -56,6 +56,68 @@ function testLoadKcatSourceQuotedRows_tc0302(testCase) verifyError(testCase, @() loadKcatSource(model, file), 'loadKcatSource:quote') end +function testReadDLKcatOutputKcatFields_tc0303(testCase) + % A kcat field that is not a finite number ('inf', '1.2.3', 'NA', 'None') + % drops the row before the reactions are checked; a file without a + % number raises. + [model, adapter, project] = kcatProject(); + file = fullfile(project, 'DLKcat_out.tsv'); + head = sprintf('rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat'); + row = @(r, s, k) sprintf('%s\tG1\t%s\tsmiles\tseq\t%s', r, s, k); + read = @() readDLKcatOutput(model, 'outFile', file, 'modelAdapter', adapter); + + writeLines(file, {head, row('R2_EXP_1', 'm1', '12.5'), row('R2_EXP_2', 'm1', 'inf'), ... + row('R2_EXP_2', 'm1', '1.2.3')}); + kcatList = read(); + verifyEqual(testCase, kcatList.rxns, {'R2_EXP_1'}) + verifyEqual(testCase, kcatList.kcats, 12.5) + + writeLines(file, {head, row('R2_EXP_1', 'm1', 'None'), row('R2_EXP_2', 'm1', 'None')}); + verifyError(testCase, read, 'readDLKcatOutput:noKcats') + + writeLines(file, {head, row('R2_EXP_1', 'm1', '12.5'), row('NOTAREACTION_EXP_1', 'm1', 'NA')}); + kcatList = read(); + verifyEqual(testCase, kcatList.rxns, {'R2_EXP_1'}) + + writeLines(file, {head, row('R2_EXP_1', 'm1', '12.5'), row('NOTAREACTION_EXP_1', 'm1', '7')}); + verifyError(testCase, read, 'readDLKcatOutput:unknownReaction') +end + +function testReadDLKcatOutputUnknownSubstrate_tc0304(testCase) + % A row with a substrate that is not in the model is left out with a + % warning; strict raises. + [model, adapter, project] = kcatProject(); + file = fullfile(project, 'DLKcat_out.tsv'); + writeLines(file, {sprintf('rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat'), ... + sprintf('R2_EXP_1\tG1\tm1\tsmiles\tseq\t12.5'), sprintf('R2_EXP_2\tG1\tNOTAREALMETABOLITE\tsmiles\tseq\t7')}); + verifyWarning(testCase, @() readDLKcatOutput(model, 'outFile', file, 'modelAdapter', adapter), ... + 'readDLKcatOutput:unknownSubstrate') + state = warning('off', 'readDLKcatOutput:unknownSubstrate'); + restore = onCleanup(@() warning(state)); + kcatList = readDLKcatOutput(model, 'outFile', file, 'modelAdapter', adapter); + verifyEqual(testCase, kcatList.rxns, {'R2_EXP_1'}) + verifyError(testCase, @() readDLKcatOutput(model, 'outFile', file, 'modelAdapter', adapter, ... + 'strict', true), 'readDLKcatOutput:unknownSubstrate') +end + +function testWriteDLKcatInputEmptySelection_tc0305(testCase) + % A selection without a reaction writes an empty file and returns no + % rows; the rows of a selection follow model.ec.rxns. + [model, adapter, project] = kcatProject(); + file = fullfile(project, 'DLKcat_in.tsv'); + none = false(numel(model.ec.rxns), 1); + evalc(['written = writeDLKcatInput(model, ''ecRxns'', none, ''modelAdapter'', adapter, ' ... + '''onlyWithSmiles'', false, ''filename'', file, ''overwrite'', true);']); + verifyEqual(testCase, size(written), [6 0]) + info = dir(file); + verifyEqual(testCase, info.bytes, 0) + some = ismember(model.ec.rxns, {'R5', 'R3'}); + evalc(['written = writeDLKcatInput(model, ''ecRxns'', some, ''modelAdapter'', adapter, ' ... + '''onlyWithSmiles'', false, ''filename'', file, ''overwrite'', true);']); + order = find(some); + verifyEqual(testCase, unique(written(1, :), 'stable'), model.ec.rxns(order)') +end + function [model, adapter, project] = kcatProject() adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); project = tempname; From 1698161804a1929c40b3bee0ec9feed1798594d2 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:12 +0100 Subject: [PATCH 19/31] kcatAssignmentConfig rejects invalid values; kcatsFromECTree takes an explicit empty ecRxns as no reaction; getStandardKcat matches pseudoreaction names in any case and runs twice on a light ecModel --- src/gather_kcats/getStandardKcat.m | 12 ++-- src/gather_kcats/kcatAssignmentConfig.m | 22 +++++-- src/gather_kcats/kcatsFromECTree.m | 21 ++++--- test/VerificationMatrix.tsv | 8 +-- test/unit_tests/kcatContractTests.m | 78 +++++++++++++++++++++++++ 5 files changed, 118 insertions(+), 23 deletions(-) diff --git a/src/gather_kcats/getStandardKcat.m b/src/gather_kcats/getStandardKcat.m index 6cad6083..d1d9c867 100644 --- a/src/gather_kcats/getStandardKcat.m +++ b/src/gather_kcats/getStandardKcat.m @@ -15,7 +15,9 @@ % first subSystem will be considered if multiple are annotated to the same % reaction. % -% Exchange, transport and pseudoreactions are filtered out, plus any +% Exchange (RAVEN's getExchangeRxns: no substrate or no product), transport +% and pseudoreactions (a reaction name that contains 'pseudoreaction' or +% 'SLIME rxn', in any letter case) are filtered out, plus any % reaction identifiers specified in /data/pseudoRxns.tsv in the model % adapter folder. % @@ -197,14 +199,14 @@ [~, exchangeRxns] = getExchangeRxns(model); transportRxns = getTransportRxns(model); [spontaneousRxns, ~] = modelAdapter.getSpontaneousReactions(model); -pseudoRxns = contains(model.rxnNames,'pseudoreaction'); -slimeRxns = contains(model.rxnNames,'SLIME rxn'); +pseudoRxns = contains(lower(model.rxnNames),'pseudoreaction'); +slimeRxns = contains(lower(model.rxnNames),'slime rxn'); rxnsToIgnore = [customRxns; exchangeRxns; find(transportRxns); ... find(spontaneousRxns); find(pseudoRxns); find(slimeRxns)]; rxnsMissingGPR(ismember(rxnsMissingGPR, rxnsToIgnore)) = []; -% Only add if not geckoLight & getStandardKcat was not run earlier -if ~any(strcmp(model.mets,'prot_standard')) +% Only add if getStandardKcat was not run earlier +if ~any(strcmp(model.mets,'prot_standard')) && ~any(strcmp(model.ec.enzymes,'standard')) % Add a new gene to be consistent with ec field named standard proteinStdGenes.genes = 'standard'; if isfield(model,'geneShortNames') diff --git a/src/gather_kcats/kcatAssignmentConfig.m b/src/gather_kcats/kcatAssignmentConfig.m index b4aca434..f6e1f57d 100644 --- a/src/gather_kcats/kcatAssignmentConfig.m +++ b/src/gather_kcats/kcatAssignmentConfig.m @@ -4,8 +4,10 @@ % Reads [kcat_assignment] and [kcat_sources.*] of model_adapter.toml (or % params.kcat_assignment and params.kcat_sources of an adapter class), and % expands a preset into its stages, as collectKcats and updateKcatSources -% use them. Without a preset or stages the preset is brenda_ec_tree. -% geckopy's kcat_assignment_config gives the same settings. +% use them. Without a preset or stages the preset is brenda_ec_tree. An +% empty preset, an unknown criteria (the names are lower case), and a stage +% whose sources is not a list raise kcatAssignment:config. geckopy's +% kcat_assignment_config gives the same settings. % % Parameters % ---------- @@ -43,7 +45,10 @@ config.isozymes = getOr(a, 'isozymes', true); config.standard = getOr(a, 'standard', true); config.sources = struct(); -hasPreset = isfield(a, 'preset') && ~isempty(a.preset); +hasPreset = isfield(a, 'preset') && ~(isnumeric(a.preset) && isempty(a.preset)); +if hasPreset && (~(ischar(a.preset) || isstring(a.preset)) || strlength(string(a.preset)) == 0) + error('kcatAssignment:config', 'The preset of [kcat_assignment] must be the name of a preset.') +end hasStages = isfield(a, 'stage') && ~isempty(a.stage); if hasPreset && hasStages error('kcatAssignment:config', ['[kcat_assignment] gives both a preset and stages; ' ... @@ -140,14 +145,19 @@ '"select" and "combine".'], stage.action) end sources = getOr(s, 'sources', {}); -if ischar(sources) - sources = {sources}; +if ~iscell(sources) + error('kcatAssignment:config', 'The sources of a [[kcat_assignment.stage]] must be a list of source names.') end if isempty(sources) error('kcatAssignment:config', 'A [[kcat_assignment.stage]] needs sources.') end stage.sources = sources(:)'; -stage.criteria = validatestring(getOr(s, 'criteria', 'max'), {'max','min','median','mean'}); +stage.criteria = getOr(s, 'criteria', 'max'); +if ~((ischar(stage.criteria) || isstring(stage.criteria)) && ... + any(strcmp(stage.criteria, {'max','min','median','mean'}))) + error('kcatAssignment:config', 'Unknown criteria of a [[kcat_assignment.stage]]; criteria: max, min, median, mean.') +end +stage.criteria = char(stage.criteria); stage.topOriginLimit = getOr(s, 'top_origin_limit', 4); stage.bottomOriginLimit = getOr(s, 'bottom_origin_limit', 6); stage.wildcardLimit = getOr(s, 'wildcard_limit', 3); diff --git a/src/gather_kcats/kcatsFromECTree.m b/src/gather_kcats/kcatsFromECTree.m index 17036801..e14828a9 100644 --- a/src/gather_kcats/kcatsFromECTree.m +++ b/src/gather_kcats/kcatsFromECTree.m @@ -21,8 +21,10 @@ % k : double % shrinkage constant (default 5): the number of measurements at which % an EC's own median and its class count equally. -% ecRxns : logical -% for which model.ec.rxns kcats are estimated (default: all). +% ecRxns : logical or double +% for which model.ec.rxns kcats are estimated, as a logical vector or as +% indices (default: all). A reaction named more than once gives one +% row. An empty ecRxns given explicitly selects no reaction. % measurements : struct % BRENDA measurement table (default: loadBrendaMeasurements with the % model adapter). @@ -59,6 +61,7 @@ 'saveAs', 'ec_tree'; ... 'modelAdapter', []}); modelAdapter = p.modelAdapter; +hasSubset = any(strcmp(cellfun(@char, varargin(1:2:end), 'UniformOutput', false), 'ecRxns')); if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); if isempty(modelAdapter) @@ -70,7 +73,7 @@ end if p.save file = fullfile(modelAdapter.params.path, 'data', 'kcats', [char(p.saveAs) '.tsv']); - warnSubset(p.ecRxns, numel(model.ec.rxns), file); + warnSubset(p.ecRxns, hasSubset, numel(model.ec.rxns), file); end measurements = p.measurements; if isempty(measurements) @@ -78,8 +81,10 @@ end idx = (1:numel(model.ec.rxns))'; -if ~isempty(p.ecRxns) - idx = idx(p.ecRxns(:)); +if hasSubset + selected = false(numel(idx), 1); + selected(p.ecRxns(:)) = true; + idx = idx(selected); end eccodes = model.ec.eccodes(idx); eccodes(cellfun(@(e) ~ischar(e) && ~isstring(e), eccodes)) = {''}; @@ -94,7 +99,7 @@ if p.save parameters = struct('k', p.k); - if ~isempty(p.ecRxns) + if hasSubset parameters.ec_rxns = countRxns(p.ecRxns); end saveKcatSource(model, kcatList, p.saveAs, 'modelAdapter', modelAdapter, ... @@ -112,9 +117,9 @@ end end -function warnSubset(ecRxns, nRxns, file) +function warnSubset(ecRxns, hasSubset, nRxns, file) % A subset written under the default name replaces the full file. -if ~isempty(ecRxns) && countRxns(ecRxns) < nRxns && isfile(file) +if hasSubset && countRxns(ecRxns) < nRxns && isfile(file) printOrange(sprintf(['%s is replaced by kcats for %d reactions only. Give saveAs ' ... 'to keep it.\n'], file, countRxns(ecRxns))); end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 22111205..6c18bb37 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -10,7 +10,7 @@ change_model/getKcatAcrossIsoenzymes.m tc0011, tc0252 change_model/getSubsetEcModel.m tc0258, tc0259 change_model/applyCustomKcats.m tc0227, tc0228, tc0229, tc0230 gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239 -gather_kcats/getStandardKcat.m tc0256, tc0257 +gather_kcats/getStandardKcat.m tc0256, tc0257, tc0309, tc0310 gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002, tc0303, tc0304 - The test doesn't really check the function, only that some output is read gather_kcats/runDLKcat.m tcm0001, tcm0002 - The test doesn't really check the function, only that we get something @@ -25,10 +25,10 @@ gather_kcats/kcatsFromOKP.m tc0075 - reads a stored result; submitting needs the gather_kcats/collectKcats.m tc0079, tc0080, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090, tc0098, tc0300 gather_kcats/assignKcats.m tc0079, tc0083 gather_kcats/updateKcatSources.m tc0082, tc0087 -gather_kcats/kcatAssignmentConfig.m tc0079, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090 +gather_kcats/kcatAssignmentConfig.m tc0079, tc0081, tc0083, tc0087, tc0088, tc0089, tc0090, tc0306 gather_kcats/loadBrendaMeasurements.m tc0085 -gather_kcats/ecTreeEstimate.m tc0084 -gather_kcats/kcatsFromECTree.m tc0086, tc0087 +gather_kcats/ecTreeEstimate.m tc0084, tc0307 +gather_kcats/kcatsFromECTree.m tc0086, tc0087, tc0308 gather_kcats/kcatsFromHomologs.m tc0095, tc0096 gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 diff --git a/test/unit_tests/kcatContractTests.m b/test/unit_tests/kcatContractTests.m index 32f19998..e61e9596 100644 --- a/test/unit_tests/kcatContractTests.m +++ b/test/unit_tests/kcatContractTests.m @@ -118,6 +118,84 @@ function testWriteDLKcatInputEmptySelection_tc0305(testCase) verifyEqual(testCase, unique(written(1, :), 'stable'), model.ec.rxns(order)') end +function testKcatAssignmentConfigRejectsInvalidValues_tc0306(testCase) + % An empty preset, a criteria that is not one of the lower case names + % and a stage whose sources is not a list raise. + config = @(a) kcatAssignmentConfig(struct('kcat_assignment', a)); + verifyError(testCase, @() config(struct('preset', '')), 'kcatAssignment:config') + verifyError(testCase, @() config(struct('stage', struct('sources', 'custom'))), 'kcatAssignment:config') + verifyError(testCase, @() config(struct('stage', struct('sources', {{'custom'}}, 'criteria', 'Median'))), ... + 'kcatAssignment:config') + verifyError(testCase, @() config(struct('stage', struct('sources', {{'custom'}}, 'criteria', 'med'))), ... + 'kcatAssignment:config') + c = config(struct('stage', struct('sources', {{'custom'}}, 'criteria', 'median'))); + verifyEqual(testCase, c.stages.criteria, 'median') + verifyEqual(testCase, config(struct()).preset, 'brenda_ec_tree') +end + +function testEcTreeEstimateRejectsNaNK_tc0307(testCase) + m = struct('ec', {{'1.1.1.1'; '1.1.1.1'}}, 'log10kcat', [1; 2]); + verifyError(testCase, @() ecTreeEstimate({'1.1.1.1'}, m, NaN), 'ecTreeEstimate:k') +end + +function testKcatsFromECTreeSubset_tc0308(testCase) + % An explicit empty ecRxns selects no reaction; a reaction named twice + % gives one row. + [model, adapter] = kcatProject(); + ecs = unique(strtrim(split(strjoin(model.ec.eccodes(~cellfun(@isempty, model.ec.eccodes)), ';'), ';'))); + ecs = ecs(~cellfun(@isempty, ecs)); + v = repelem((1:numel(ecs))', 2); + m = struct('ec', {repelem(ecs, 2)}, 'log10kcat', v, 'kcat', 10 .^ v, 'release', 'BRENDA release test'); + withEc = find(~cellfun(@isempty, model.ec.eccodes)); + verifyGreaterThanOrEqual(testCase, numel(withEc), 2) + everything = kcatsFromECTree(model, 'measurements', m, 'save', false, 'modelAdapter', adapter); + verifyEqual(testCase, numel(everything.rxns), numel(withEc)) + none = kcatsFromECTree(model, 'measurements', m, 'save', false, 'modelAdapter', adapter, 'ecRxns', []); + verifyEmpty(testCase, none.rxns) + twice = kcatsFromECTree(model, 'measurements', m, 'save', false, 'modelAdapter', adapter, ... + 'ecRxns', [withEc(1); withEc(2); withEc(2)]); + verifyEqual(testCase, twice.rxns, model.ec.rxns(withEc(1:2))) +end + +function testGetStandardKcatExchangeAndPseudoReactions_tc0309(testCase) + % R7 consumes two metabolites and produces none: an exchange reaction, + % without the standard enzyme. R8 is named like a pseudoreaction in + % another letter case. R6 is an ordinary reaction without a GPR. + adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); + model = loadConventionalGEM('modelAdapter', adapter); + toAdd.rxns = {'R6'; 'R7'; 'R8'}; + toAdd.rxnNames = {'R6'; 'R7'; 'Biomass Pseudoreaction'}; + toAdd.equations = {'m1[c] => m2[c]'; 'm1[c] + m2[c] =>'; 'm1[c] => m2[c]'}; + toAdd.grRules = {''; ''; ''}; + model = addRxns(model, toAdd, 3); + ecModel = makeEcModel(model, false, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 10; + evalc('ecModel = getStandardKcat(ecModel, ''modelAdapter'', adapter);'); + verifyTrue(testCase, ismember('R6', ecModel.ec.rxns)) + verifyFalse(testCase, ismember('R7', ecModel.ec.rxns)) + verifyFalse(testCase, ismember('R8', ecModel.ec.rxns)) +end + +function testGetStandardKcatLightModelTwice_tc0310(testCase) + % The standard pseudoenzyme of a light model is not added again, and the + % second call gives the first call's result. + adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); + model = loadConventionalGEM('modelAdapter', adapter); + toAdd.rxns = {'R6'}; + toAdd.equations = {'m1[c] => m2[c]'}; + toAdd.grRules = {''}; + model = addRxns(model, toAdd, 3); + ecModel = makeEcModel(model, true, adapter); + ecModel = getECfromGEM(ecModel); + ecModel.ec.kcat(:) = 10; + evalc('first = getStandardKcat(ecModel, ''modelAdapter'', adapter);'); + evalc('second = getStandardKcat(first, ''modelAdapter'', adapter);'); + verifyEqual(testCase, sum(strcmp(second.ec.enzymes, 'standard')), 1) + verifyEqual(testCase, sort(second.ec.rxns), sort(first.ec.rxns)) + verifyEqual(testCase, sort(second.ec.kcat), sort(first.ec.kcat)) +end + function [model, adapter, project] = kcatProject() adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); project = tempname; From 32a49d725f2625f7ddf7ded60d42de25f04658e4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 21:55:33 +0100 Subject: [PATCH 20/31] saveEcModel writes a default metaData description when there is none A model without annotation.description gets 'Enzyme-constrained model of '; an existing description is kept, as save_ec_model does. --- src/io/saveEcModel.m | 6 ++++++ test/VerificationMatrix.tsv | 2 +- test/unit_tests/parityRound5G4Tests.m | 19 +++++++++++++++++++ 3 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/io/saveEcModel.m b/src/io/saveEcModel.m index 9b1d8f65..26116d08 100644 --- a/src/io/saveEcModel.m +++ b/src/io/saveEcModel.m @@ -10,6 +10,9 @@ function saveEcModel(ecModel,varargin) % A YAML file records its ecModel format (gecko_format, see loadEcModel) % and the GECKO release that wrote it (generator) in its metaData, in % place of the geckopy_version that earlier geckopy releases wrote. +% A model without a description (ecModel.annotation.description) gets +% 'Enzyme-constrained model of ' in its metaData; an existing description +% is kept. % % The ecModel is saved to the ecModel-specific models/ subfolder. For % saving to other locations, writeYAMLmodel or exportModel are more @@ -88,6 +91,9 @@ function saveEcModel(ecModel,varargin) % model_adapter.toml) and the release that wrote it. ecModel.annotation.gecko_format = '4'; ecModel.annotation.generator = ['GECKO ' geckoVersion()]; + if ~isfield(ecModel.annotation, 'description') || isempty(ecModel.annotation.description) + ecModel.annotation.description = ['Enzyme-constrained model of ' ecModel.id]; + end if isfield(ecModel.annotation, 'geckopy_version') ecModel.annotation = rmfield(ecModel.annotation, 'geckopy_version'); end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 07346755..389c9a49 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -79,7 +79,7 @@ model_adapter/adapterTemplate.toml tc0068 utilities/* Not explicitly tested - The code has existed for a long time. There are a lot of functions in here, unclear which of them are still useful and functional. io/loadConventionalGEM.m tc0009 io/loadEcModel.m tc0009, tc0099 -io/saveEcModel.m tc0009, tc0099 +io/saveEcModel.m tc0009, tc0099, tc0603 tutorials/tutorial_Human-GEM/HumanGEMAdapter.m tcm0002 - does not test all aspects of the adapter, but some. tutorials/full_ecModel/model_adapter.toml tcm0001 - does not test all aspects of the adapter, but some. tutorials/tutorial_yeast-GEM/code/changeMedia_batch.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/parityRound5G4Tests.m b/test/unit_tests/parityRound5G4Tests.m index 598a27b2..e6d07c28 100644 --- a/test/unit_tests/parityRound5G4Tests.m +++ b/test/unit_tests/parityRound5G4Tests.m @@ -33,6 +33,25 @@ function testFillCarbonNumUsesFormulaCarbons_tc0602(testCase) verifyEqual(testCase, out.excarbon(strcmp(out.rxns,'EX_a')), 6) end +function testSaveEcModelWritesDefaultDescriptionOnlyWhenNone_tc0603(testCase) + adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); + project = tempname; + testCase.addTeardown(@() rmdir(project, 's')); + copyfile(fullfile(adapter.params.path, 'data'), fullfile(project, 'data')); + adapter.params.path = project; + mkdir(fullfile(project, 'models')); + model = getECfromGEM(makeEcModel(getGeckoTestModel(), false, adapter)); + model.annotation.description = ''; + saveEcModel(model, 'default.yml', adapter); + text = fileread(fullfile(project, 'models', 'default.yml')); + verifyNotEmpty(testCase, strfind(text, ['- description: Enzyme-constrained model of ' model.id])) + model.annotation.description = 'my own description'; + saveEcModel(model, 'own.yml', adapter); + text = fileread(fullfile(project, 'models', 'own.yml')); + verifyNotEmpty(testCase, strfind(text, '- description: my own description')) + verifyEmpty(testCase, strfind(text, 'Enzyme-constrained model of')) +end + function model = exchangeToyModel() model = struct(); model.id = 'toy'; From 97e88941d5d2cb755f10977b037d3b3589b0a0c0 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 21:57:51 +0100 Subject: [PATCH 21/31] applyKcatConstraints takes cobra and ec reaction ids on a light model and clears its pool coefficients when every kcat is 0 --- src/change_model/applyKcatConstraints.m | 39 ++++++++++++++++++------- test/VerificationMatrix.tsv | 2 +- test/unit_tests/parityRound3Tests.m | 13 +++++---- test/unit_tests/parityRound5G3Tests.m | 22 ++++++++++++++ 4 files changed, 59 insertions(+), 17 deletions(-) diff --git a/src/change_model/applyKcatConstraints.m b/src/change_model/applyKcatConstraints.m index e87f7ce2..1edabd5e 100644 --- a/src/change_model/applyKcatConstraints.m +++ b/src/change_model/applyKcatConstraints.m @@ -3,9 +3,8 @@ % % Existing enzyme constraints are first removed (unless updateRxns is % provided), and new constraints are defined based on the content of -% model.ec.kcat. When every kcat is 0, the existing constraints of a full -% ecModel are removed and no new ones are written; a light ecModel is -% returned unchanged. +% model.ec.kcat. When every kcat is 0, the existing constraints are removed +% and no new ones are written. % % Parameters % ---------- @@ -18,8 +17,13 @@ % if not all enzyme constraints should be updated, this can be given as % either a logical vector of length model.ec.rxns, a vector of % model.ec.rxns indices, or a (cell array of) string(s) with -% model.ec.rxns identifiers. For light models, these reactions should -% match model.rxns. +% model.ec.rxns identifiers. For light models, a (cell array of) string(s) +% may hold reaction identifiers of model.rxns or of model.ec.rxns (with +% the ###_ prefix); a reaction is then updated as a whole, with all its +% isozymes. An identifier that is in both or in neither raises the +% error applyKcatConstraints:ambiguousReaction or +% applyKcatConstraints:unknownReaction. A logical vector or indices +% refer to model.rxns for light models. % % Returns % ------- @@ -56,6 +60,23 @@ updateRxns = updateRxnsLog; elseif iscellstr(updateRxns) || ischar(updateRxns) || isstring(updateRxns) updateRxnsIds = convertCharArray(updateRxns); + updateRxnsIds = updateRxnsIds(:); + if model.ec.geckoLight + isCobra = ismember(updateRxnsIds,model.rxns); + isEc = ismember(updateRxnsIds,model.ec.rxns); + if any(isCobra & isEc) + error('applyKcatConstraints:ambiguousReaction', ... + 'The identifiers %s are in both model.rxns and model.ec.rxns.', ... + strjoin(updateRxnsIds(isCobra & isEc),', ')) + end + if any(~isCobra & ~isEc) + error('applyKcatConstraints:unknownReaction', ... + 'The identifiers %s are in neither model.rxns nor model.ec.rxns.', ... + strjoin(updateRxnsIds(~isCobra & ~isEc),', ')) + end + updateRxnsIds = [updateRxnsIds(isCobra); extractAfter(updateRxnsIds(isEc),4)]; + updateRxnsIds = updateRxnsIds(:); + end updateRxns = ismember(rxns,updateRxnsIds); end @@ -68,12 +89,8 @@ ' not a valid GECKO3 ecModel. First run makeEcModel(model).']) end if all(model.ec.kcat==0) - if model.ec.geckoLight - printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n'); - return - end - % In a full model the existing enzyme constraints of the selected - % reactions are still cleared below, and no new ones are written. + % The existing enzyme constraints of the selected reactions are still + % cleared below, and no new ones are written. printOrange('WARNING: No kcat values are provided in model.ec.kcat, existing enzyme constraints are cleared and no new ones are written.\n'); end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 5136f65f..07dfc347 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -1,7 +1,7 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time change_model/applyComplexData.m tc0003, tc0004, tc0103, tc0501 -change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238, tc0250 +change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238, tc0250, tc0512 change_model/findMetSmiles.m tc0012 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 diff --git a/test/unit_tests/parityRound3Tests.m b/test/unit_tests/parityRound3Tests.m index b28b74e8..fa9edb33 100644 --- a/test/unit_tests/parityRound3Tests.m +++ b/test/unit_tests/parityRound3Tests.m @@ -22,9 +22,9 @@ function testApplyKcatConstraintsAllZeroClearsFullModel_tc0250(testCase) verifyEqual(testCase, cleared.S(:, ~metabolic), ecModel.S(:, ~metabolic)) end -function testApplyKcatConstraintsAllZeroLightIsUnchanged_tc0250light(testCase) - % A light ecModel with every kcat 0 is returned unchanged (control: the - % behaviour is the same before and after the full-model change). +function testApplyKcatConstraintsAllZeroClearsLightModel_tc0250light(testCase) + % With every kcat 0 the pool coefficients of a light ecModel are cleared, + % as the enzyme constraints of a full ecModel are. adapter = getAdapter(); lecModel = makeEcModel(getGeckoTestModel(), true, adapter); lecModel.ec.kcat(:) = 10; @@ -32,8 +32,11 @@ function testApplyKcatConstraintsAllZeroLightIsUnchanged_tc0250light(testCase) poolRow = strcmp(lecModel.mets, 'prot_pool'); verifyGreaterThan(testCase, nnz(lecModel.S(poolRow, :)), 1) lecModel.ec.kcat(:) = 0; - same = applyKcatConstraints(lecModel); - verifyEqual(testCase, same.S, lecModel.S) + cleared = applyKcatConstraints(lecModel); + poolRow = strcmp(cleared.mets, 'prot_pool'); + metabolic = ~strcmp(cleared.rxns, 'prot_pool_exchange'); + verifyEqual(testCase, nnz(cleared.S(poolRow, metabolic)), 0) + verifyEqual(testCase, cleared.S(:, ~metabolic), lecModel.S(:, ~metabolic)) end function testSensitivityTuningStopsOnGrowthBelowTolerance_tc0251(testCase) diff --git a/test/unit_tests/parityRound5G3Tests.m b/test/unit_tests/parityRound5G3Tests.m index a9097005..a37fd604 100644 --- a/test/unit_tests/parityRound5G3Tests.m +++ b/test/unit_tests/parityRound5G3Tests.m @@ -208,6 +208,28 @@ function testRelaxProteomicsGreedyEmptyEnzymeSet_tc0511(testCase) verifyEqual(testCase, fieldnames(all.relaxed), {'P5'}) end +function testApplyKcatConstraintsLightUpdateIds_tc0512(testCase) + % On a light model the update list takes cobra reaction ids and ec + % reaction ids (with the ###_ prefix); an id of neither kind raises. + adapter = getAdapter(); + lecModel = makeEcModel(getGeckoTestModel(), true, adapter); + lecModel.ec.kcat(:) = 10; + poolRow = strcmp(lecModel.mets, 'prot_pool'); + all = applyKcatConstraints(lecModel); + stale = all; + stale.S(poolRow, strcmp(stale.rxns, 'R3')) = 0; + byCobra = applyKcatConstraints(stale, 'updateRxns', {'R3'}); + verifyEqual(testCase, byCobra.S, all.S) + ecId = lecModel.ec.rxns(endsWith(lecModel.ec.rxns, '_R3')); + byEc = applyKcatConstraints(stale, 'updateRxns', ecId); + verifyEqual(testCase, byEc.S, all.S) + verifyError(testCase, @() applyKcatConstraints(lecModel, 'updateRxns', {'NOT_A_RXN'}), ... + 'applyKcatConstraints:unknownReaction') + lecModel.rxns{end+1} = lecModel.ec.rxns{1}; + verifyError(testCase, @() applyKcatConstraints(lecModel, 'updateRxns', lecModel.ec.rxns(1)), ... + 'applyKcatConstraints:ambiguousReaction') +end + function adapter = getAdapter() geckoPath = findGECKOroot; adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); From 73c52b45b8cd12f679d9e03a8ed2aadb67c36798 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 21:58:40 +0100 Subject: [PATCH 22/31] assignKcatValues writes ec.source as the lower case source followed by the wildcard level and origin in brackets, as geckopy does --- src/gather_kcats/assignKcatValues.m | 33 ++++++++++++++++++++++++++- test/VerificationMatrix.tsv | 2 +- test/unit_tests/kcatAssignmentTests.m | 2 +- test/unit_tests/kcatContractTests.m | 13 +++++++++++ 4 files changed, 47 insertions(+), 3 deletions(-) diff --git a/src/gather_kcats/assignKcatValues.m b/src/gather_kcats/assignKcatValues.m index 8188a3b0..227ca798 100644 --- a/src/gather_kcats/assignKcatValues.m +++ b/src/gather_kcats/assignKcatValues.m @@ -33,6 +33,11 @@ % 'median'/'mean', where the selected value need not equal any single % input, it is taken from the first entry instead (arbitrary but % deterministic, since an aggregate has no single "winning" source). +% The source is written in lower case snake_case, with the wildcard +% level and origin of the entry in brackets when it has them +% (kcatList.wildcardLvl and kcatList.origin), for example +% 'brenda (wc=0, origin=1)'. Functions that compare model.ec.source with +% a source name use the part before ' (' (see resolveTunableContext). % overwrite : char % whether existing kcat values should be overwritten. Options: 'true', % 'false', 'ifHigher'. The last option will overwrite only if the new @@ -77,6 +82,17 @@ % Remove zero kcat values. Only adjusting fields that are used later. removeZero = kcatList.kcats == 0; +nEntries = numel(removeZero); +wildcardLvl = nan(nEntries, 1); +origin = nan(nEntries, 1); +if isfield(kcatList,'wildcardLvl') && numel(kcatList.wildcardLvl) == nEntries + wildcardLvl = double(kcatList.wildcardLvl(:)); +end +if isfield(kcatList,'origin') && numel(kcatList.origin) == nEntries + origin = double(kcatList.origin(:)); +end +wildcardLvl(removeZero) = []; +origin(removeZero) = []; kcatList.kcats(removeZero) = []; kcatList.rxns(removeZero) = []; if isfield(kcatList,'kcatSource') @@ -114,7 +130,8 @@ otherwise error('Invalid criteria specified') end - selectedSource(i) = kcatList.kcatSource(idxMatch(j)); + selectedSource{i} = formatKcatSource(kcatList.kcatSource{idxMatch(j)}, ... + wildcardLvl(idxMatch(j)), origin(idxMatch(j))); end % Populate model.ec.kcat @@ -140,3 +157,17 @@ end rxnIdx = idxInModelUnique; end + +function text = formatKcatSource(source, wildcardLvl, origin) +% The provenance text of model.ec.source: the source label in lower case +% snake_case, followed by the wildcard level and origin of a fuzzy BRENDA +% match in brackets, e.g. 'brenda (wc=0, origin=1)'. +text = lower(strtrim(char(source))); +text = regexprep(text, '[^0-9a-z]+', '_'); +text = regexprep(text, '^_+|_+$', ''); +if ~isnan(wildcardLvl) && ~isnan(origin) + text = sprintf('%s (wc=%d, origin=%d)', text, fix(wildcardLvl), fix(origin)); +elseif ~isnan(wildcardLvl) + text = sprintf('%s (wc=%d)', text, fix(wildcardLvl)); +end +end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 6c18bb37..2d50952a 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -14,7 +14,7 @@ gather_kcats/getStandardKcat.m tc0256, tc0257, tc0309, tc0310 gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241 gather_kcats/readDLKcatOutput.m tcm0001, tcm0002, tc0303, tc0304 - The test doesn't really check the function, only that some output is read gather_kcats/runDLKcat.m tcm0001, tcm0002 - The test doesn't really check the function, only that we get something -gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameterizations here; also tc0104 +gather_kcats/assignKcatValues.m tc0011 - we did not test all possible parameterizations here; also tc0104, tc0311 gather_kcats/writeDLKcatInput.m partly by tc0011, tc0305 gather_kcats/getEcRxnKeys.m tc0070, tc0076 gather_kcats/saveKcatSource.m tc0071, tc0073, tc0076, tc0078 diff --git a/test/unit_tests/kcatAssignmentTests.m b/test/unit_tests/kcatAssignmentTests.m index 68312c88..97cdae66 100644 --- a/test/unit_tests/kcatAssignmentTests.m +++ b/test/unit_tests/kcatAssignmentTests.m @@ -63,7 +63,7 @@ function testStagesFromToml_tc0081(testCase) adapter.params.kcat_assignment = struct('stage', struct('action', 'select', ... 'sources', {{'dlkcat'}}, 'criteria', 'min'), 'isozymes', false, 'standard', false); decisions = collectKcats(base, 'modelAdapter', adapter, 'writeReport', false); - verifyEqual(testCase, unique(decisions.source(decisions.kcat > 0)), {'DLKcat'}) + verifyEqual(testCase, unique(decisions.source(decisions.kcat > 0)), {'dlkcat'}) verifyEqual(testCase, decisions.kcat(strcmp(decisions.ecRxn, 'R2_EXP_1')), 6.5) adapter.params.kcat_assignment.preset = 'gecko_median'; verifyError(testCase, @() collectKcats(base, 'modelAdapter', adapter), 'kcatAssignment:config') diff --git a/test/unit_tests/kcatContractTests.m b/test/unit_tests/kcatContractTests.m index e61e9596..d6a0db87 100644 --- a/test/unit_tests/kcatContractTests.m +++ b/test/unit_tests/kcatContractTests.m @@ -196,6 +196,19 @@ function testGetStandardKcatLightModelTwice_tc0310(testCase) verifyEqual(testCase, sort(second.ec.kcat), sort(first.ec.kcat)) end +function testAssignKcatValuesWritesSourceWithMatchDetail_tc0311(testCase) + % ec.source is the source in lower case snake_case, followed by the + % wildcard level and origin of a BRENDA match in brackets, as geckopy + % writes it. + [model, ~] = kcatProject(); + rxns = model.ec.rxns(1:4); + kcatList = struct('source', 'x', 'rxns', {rxns}, 'kcats', [5; 6; 7; 8], ... + 'kcatSource', {{'BRENDA'; 'DLKcat'; 'OKP-CataPro'; 'brenda'}}, ... + 'wildcardLvl', [0; NaN; NaN; 2], 'origin', [1; NaN; 3; NaN]); + model = assignKcatValues(model, kcatList); + verifyEqual(testCase, model.ec.source(1:4), {'brenda (wc=0, origin=1)'; 'dlkcat'; 'okp_catapro'; 'brenda (wc=2)'}) +end + function [model, adapter, project] = kcatProject() adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); project = tempname; From 482e4cd47ebd527d5c3ce438ebf3fb84adce9a6d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 21:58:40 +0100 Subject: [PATCH 23/31] resolveTunableContext puts a kcat in the group of its source name before the bracketed detail --- src/kcat_tuning/evotune/resolveTunableContext.m | 11 ++++++++++- test/VerificationMatrix.tsv | 2 +- test/unit_tests/kcatContractTests.m | 16 ++++++++++++++++ 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/src/kcat_tuning/evotune/resolveTunableContext.m b/src/kcat_tuning/evotune/resolveTunableContext.m index 28894f14..abddb708 100644 --- a/src/kcat_tuning/evotune/resolveTunableContext.m +++ b/src/kcat_tuning/evotune/resolveTunableContext.m @@ -5,6 +5,12 @@ % Shared by screenKcatLeverage and cmaesKcatTuning so a screen and the % tuning run it feeds are always looking at exactly the same parameters. % +% A kcat belongs to the group whose params.kcatSources entry equals its +% ec.source, or the part of ec.source before ' (' (case-insensitive), so +% 'brenda (wc=0, origin=1)' is in the group of 'brenda' and +% 'combined (sd=0.10, n=2)' in the group of 'combined', with that group's +% sigma0logSource. +% % Parameters % ---------- % ecModel : struct @@ -48,11 +54,14 @@ ecRxnIdsTunable = ecModel.ec.rxns(tunableIdx); kcat0 = ecModel.ec.kcat(tunableIdx); sources = ecModel.ec.source(tunableIdx); +% A source with a detail, such as 'combined (sd=0.10, n=2)', is also matched by +% its name before ' ('. +sourceNames = strtrim(regexprep(sources, ' \(.*$', '')); groups = repmat({'unlabelled'}, numel(tunableIdx), 1); sigma0log = params.sigma0logDefault * ones(numel(tunableIdx),1); for i = 1:numel(params.kcatSources) - idx = strcmpi(sources, params.kcatSources{i}); + idx = strcmpi(sources, params.kcatSources{i}) | strcmpi(sourceNames, params.kcatSources{i}); groups(idx) = params.kcatSources(i); sigma0log(idx) = params.sigma0logSource(i); end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 2d50952a..997be86d 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -51,7 +51,7 @@ get_enzyme_data/getECstring.m Called from findECInDB - tested when that function download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 download_databases/fetchBrendaAsset.m tc0085 download_databases/loadDatabases.m Called from getECfromDatabase - tested when that function is tested. -kcat_tuning/evotune/* Has only been tested by running it manually, not tested by test cases +kcat_tuning/evotune/* resolveTunableContext by tc0312; the rest has only been tested by running it manually, not tested by test cases kcat_tuning/fillCarbonNum.m tc0055 (deprecated-alias check only, not a real usage test) kcat_tuning/findTopLimitations.m Not explicitly tested - The code has existed for a long time kcat_tuning/sensitivityTuning.m tcm0001, tc0251, tc0253, tc0254 diff --git a/test/unit_tests/kcatContractTests.m b/test/unit_tests/kcatContractTests.m index d6a0db87..1539fe32 100644 --- a/test/unit_tests/kcatContractTests.m +++ b/test/unit_tests/kcatContractTests.m @@ -209,6 +209,22 @@ function testAssignKcatValuesWritesSourceWithMatchDetail_tc0311(testCase) verifyEqual(testCase, model.ec.source(1:4), {'brenda (wc=0, origin=1)'; 'dlkcat'; 'okp_catapro'; 'brenda (wc=2)'}) end +function testResolveTunableContextMatchesSourceBeforeDetail_tc0312(testCase) + % A source with a detail is in the group of its name: 'brenda (wc=0, + % origin=1)' in brenda and 'combined (sd=0.10, n=2)' in combined, with + % that group's width; a source that no group names keeps the default. + [model, ~] = kcatProject(); + model.ec.kcat(1:4) = 10; + model.ec.source(1:4) = {'brenda (wc=0, origin=1)'; 'combined (sd=0.10, n=2)'; 'DLKcat'; 'other (x)'}; + params = struct('kcatSources', {{'brenda', 'combined', 'dlkcat'}}, 'sigma0logSource', [0.2; 0.3; 0.4], ... + 'sigma0logDefault', 0.5, 'tieIsozymes', false); + mask = false(numel(model.ec.rxns), 1); + mask(1:4) = true; + [~, ~, ~, groups, sigma0log] = resolveTunableContext(model, params, mask); + verifyEqual(testCase, groups, {'brenda'; 'combined'; 'dlkcat'; 'unlabelled'}) + verifyEqual(testCase, sigma0log, [0.2; 0.3; 0.4; 0.5]) +end + function [model, adapter, project] = kcatProject() adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot, 'test', 'unit_tests', 'ecTestGEM', 'TestGEMAdapter.m')); project = tempname; From 1903908217896c71dfec292fb7a04cc440774840 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:00:55 +0100 Subject: [PATCH 24/31] findMetSmiles uses the last row of a duplicated cache name and skips unnamed metabolites with a warning --- src/change_model/findMetSmiles.m | 14 +++++++--- test/unit_tests/metSmilesCacheTests.m | 39 +++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 3 deletions(-) create mode 100644 test/unit_tests/metSmilesCacheTests.m diff --git a/src/change_model/findMetSmiles.m b/src/change_model/findMetSmiles.m index 616f1caf..55f4f94b 100644 --- a/src/change_model/findMetSmiles.m +++ b/src/change_model/findMetSmiles.m @@ -3,7 +3,10 @@ % % Matches are cached in tutorials/***/data/smilesDB.tsv, which is checked % first on subsequent runs so repeated queries are avoided. If the model -% already has a metSmiles field, non-empty entries are left unchanged. +% already has a metSmiles field, non-empty entries are left unchanged. When +% the cache holds a name more than once, the last row is used. A metabolite +% without a name is skipped with a warning (findMetSmiles:emptyName) and +% gets no SMILES. % % Parameters % ---------- @@ -39,6 +42,10 @@ end params = modelAdapter.params; +nEmpty = sum(cellfun(@isempty,model.metNames)); +if nEmpty > 0 + warning('findMetSmiles:emptyName','%d metabolite(s) have no name and get no SMILES.',nEmpty) +end [uniqueNames, ~, uniqueIdx] = unique(regexprep(model.metNames,'^prot_.*','')); uniqueSmiles(1:numel(uniqueNames),1) = {''}; metMatch = false(length(uniqueNames),1); @@ -49,8 +56,9 @@ fID = fopen(smilesDBfile,'r'); raw = textscan(fID,'%s %s','Delimiter','\t','HeaderLines',0); fclose(fID); - smilesDB.names = raw{1}; - smilesDB.smile = raw{2}; + %The last row of a name that occurs more than once is used + smilesDB.names = flipud(raw{1}); + smilesDB.smile = flipud(raw{2}); [metMatch, metIdx] = ismember(uniqueNames,smilesDB.names); uniqueSmiles(metMatch) = smilesDB.smile(metIdx(metMatch)); if verbose; fprintf('done.\n'); end diff --git a/test/unit_tests/metSmilesCacheTests.m b/test/unit_tests/metSmilesCacheTests.m new file mode 100644 index 00000000..97391884 --- /dev/null +++ b/test/unit_tests/metSmilesCacheTests.m @@ -0,0 +1,39 @@ +%run this test case with the command +%results = runtests('metSmilesCacheTests.m') + +function tests = metSmilesCacheTests + tests = functiontests(localfunctions); +end + +function testFindMetSmilesUsesLastRowOfDuplicateName_tc0407(testCase) + adapter = cacheAdapter(testCase, sprintf('a\tFIRST\nb\tCCO\na\tLAST\n')); + model = struct('metNames', {{'a';'b';'a'}}); + [~, out] = evalc('findMetSmiles(model, ''modelAdapter'', adapter, ''verbose'', false)'); + verifyEqual(testCase, out.metSmiles, {'LAST';'CCO';'LAST'}) +end + +function testFindMetSmilesSkipsEmptyNameWithWarning_tc0408(testCase) + adapter = cacheAdapter(testCase, sprintf('a\tCCO\n')); + model = struct('metNames', {{'a';'';'prot_P1'}}); + verifyWarning(testCase, @() quietFind(model, adapter), 'findMetSmiles:emptyName') + warning('off', 'findMetSmiles:emptyName') + cleanup = onCleanup(@() warning('on', 'findMetSmiles:emptyName')); %#ok + [~, out] = evalc('findMetSmiles(model, ''modelAdapter'', adapter, ''verbose'', false)'); + verifyEqual(testCase, out.metSmiles, {'CCO';'';''}) +end + +function quietFind(model, adapter) + evalc('findMetSmiles(model, ''modelAdapter'', adapter, ''verbose'', false)'); +end + +function adapter = cacheAdapter(testCase, cacheText) + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); + folder = tempname(); + mkdir(fullfile(folder, 'data')); + testCase.addTeardown(@() rmdir(folder, 's')) + fid = fopen(fullfile(folder, 'data', 'smilesDB.tsv'), 'w'); + fprintf(fid, '%s', cacheText); + fclose(fid); + adapter.params.path = folder; +end From 4ea4faeb1d4f0ebd4acec2825642d630b12422e6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:00:56 +0100 Subject: [PATCH 25/31] getECfromDatabase raises for actions other than display and ignore --- src/get_enzyme_data/getECfromDatabase.m | 9 +++++++-- test/unit_tests/ecFromDatabaseConflictTests.m | 6 ++++++ 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/get_enzyme_data/getECfromDatabase.m b/src/get_enzyme_data/getECfromDatabase.m index 102fe838..ce498249 100644 --- a/src/get_enzyme_data/getECfromDatabase.m +++ b/src/get_enzyme_data/getECfromDatabase.m @@ -28,8 +28,8 @@ % uses the protein with the lowest index in the database. % - 'ignore' ignores multiplicities and uses the protein with the % lowest index in the database. -% - 'add' adds all the multiple proteins as isozymes for the given -% reaction. +% +% Any other action raises getECfromDatabase:unsupportedAction. % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % @@ -56,6 +56,11 @@ action = p.action; modelAdapter = p.modelAdapter; +if ~(ischar(action) && any(strcmpi(action,{'display','ignore'}))) + error('getECfromDatabase:unsupportedAction', ... + 'action must be ''display'' or ''ignore''.') +end + if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); if isempty(modelAdapter) diff --git a/test/unit_tests/ecFromDatabaseConflictTests.m b/test/unit_tests/ecFromDatabaseConflictTests.m index a511106f..1eb6f414 100644 --- a/test/unit_tests/ecFromDatabaseConflictTests.m +++ b/test/unit_tests/ecFromDatabaseConflictTests.m @@ -41,6 +41,12 @@ function testGetECfromDatabaseEmptySelectionChangesNothing_tc0406(testCase) verifyEqual(testCase, all.ec.eccodes{1}, '1.1.1.1') end +function testGetECfromDatabaseRejectsUnsupportedAction_tc0409(testCase) + [ecModel, adapter] = conflictModel(testCase); + verifyError(testCase, @() getECfromDatabase(ecModel, [], 'add', adapter), 'getECfromDatabase:unsupportedAction') + verifyError(testCase, @() getECfromDatabase(ecModel, [], 'other', adapter), 'getECfromDatabase:unsupportedAction') +end + function [ecModel, adapter] = conflictModel(testCase) % ecTestGEM with a second UniProt entry for gene G4 and another EC number. geckoPath = findGECKOroot; From 776b79fa2ce0eb62bb006f64a592c38aadd7e840 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:00:56 +0100 Subject: [PATCH 26/31] reviewKcatAssignment help text states the case rule for sources --- src/gather_kcats/reviewKcatAssignment.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gather_kcats/reviewKcatAssignment.m b/src/gather_kcats/reviewKcatAssignment.m index d53cf600..a5d265af 100644 --- a/src/gather_kcats/reviewKcatAssignment.m +++ b/src/gather_kcats/reviewKcatAssignment.m @@ -26,7 +26,8 @@ % ---------- % rxnIds, kcats, sources, ecCodes : cell/double % ec.rxns, ec.kcat, ec.source, ec.eccodes (or the equivalent -% subset), all the same length. +% subset), all the same length. The source labels 'custom' and +% 'standard' are recognised without regard to letter case. % % Name-Value Arguments % -------------------- From 5437f4ab2aea1edc3684dc0e714b7873c71cc631 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:00:56 +0100 Subject: [PATCH 27/31] Record tc0407 to tc0409 in the verification matrix --- test/VerificationMatrix.tsv | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 3244566f..10fbc6e5 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -2,7 +2,7 @@ Function/File (in src) Verification Comments BrendaParser/* Not explicitly tested - The code has existed for a long time change_model/applyComplexData.m tc0003, tc0004, tc0103 change_model/applyKcatConstraints.m tc0011, tc0215, tc0235, tc0238, tc0250 -change_model/findMetSmiles.m tc0012 +change_model/findMetSmiles.m tc0012, tc0407, tc0408 change_model/getComplexData.m tc0003, tc0004 - File download not covered by test cases but is tested manually and works. change_model/makeEcModel.m tc0001, tc0002, tc0210, tc0212, tc0233, tc0234 change_model/setProtPoolSize.m tc0005 @@ -45,7 +45,7 @@ gather_kcats/ecStatsFromBrenda.m tc0400, tc0401 - not verified against a real BR gather_kcats/mergeKcats.m tc0240 get_enzyme_data/calculateMW.m Not explicitly tested - The code has existed for a long time get_enzyme_data/findECInDB.m Called from getECfromDatabase - tested when that function is tested; also tc0242, tc0243, tc0244, tc0245, tc0246 -get_enzyme_data/getECfromDatabase.m tc0007, tc0404, tc0405, tc0406 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM +get_enzyme_data/getECfromDatabase.m tc0007, tc0404, tc0405, tc0406, tc0409 - does not test download of the databases - this is tested in the manual workflows for Yeast-GEM and Human-GEM get_enzyme_data/getECfromGEM.m tc0006 get_enzyme_data/getECstring.m Called from findECInDB - tested when that function is tested; also tc0246 download_databases/loadBRENDAdata.m Called from fuzzyKcatMatching - tested when that function is tested; also tc0105, tc0220, tc0221, tc0222 From e0b2b1e8c7fb46940726aa735f6ac3222e6563a6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:02:31 +0100 Subject: [PATCH 28/31] tc0009: the reloaded annotation carries the default description --- test/unit_tests/geckoCoreFunctionTests.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 308b830f..a6c227c5 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -242,7 +242,7 @@ function testsaveECModel_tc0009(testCase) delete(fullfile(adapter.params.path,'models','ecModel.yml')); verifyEqual(testCase, loadedEcModel.annotation.gecko_format, '4') verifyTrue(testCase, startsWith(loadedEcModel.annotation.generator, 'GECKO ')) - loadedEcModel.annotation = rmfield(loadedEcModel.annotation, {'gecko_format', 'generator'}); + loadedEcModel.annotation = rmfield(loadedEcModel.annotation, {'gecko_format', 'generator', 'description'}); verifyEqual(testCase, ecModel, loadedEcModel) % Test loading of conventional GEM From eb727bd1f057645fa65a8614f41c5a227c96c8a8 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:04:38 +0100 Subject: [PATCH 29/31] saveEcModel: a note counts as the description; tc0009 expects the note --- src/io/saveEcModel.m | 13 +++++++++---- test/unit_tests/geckoCoreFunctionTests.m | 2 +- test/unit_tests/parityRound5G4Tests.m | 6 ++++++ 3 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/io/saveEcModel.m b/src/io/saveEcModel.m index 26116d08..5d96cbc4 100644 --- a/src/io/saveEcModel.m +++ b/src/io/saveEcModel.m @@ -10,9 +10,10 @@ function saveEcModel(ecModel,varargin) % A YAML file records its ecModel format (gecko_format, see loadEcModel) % and the GECKO release that wrote it (generator) in its metaData, in % place of the geckopy_version that earlier geckopy releases wrote. -% A model without a description (ecModel.annotation.description) gets -% 'Enzyme-constrained model of ' in its metaData; an existing description -% is kept. +% A model without a description gets 'Enzyme-constrained model of ' in +% its metaData; an existing description is kept. readYAMLmodel reads the +% metaData description into ecModel.annotation.note, so a note counts as the +% description of the model, as does ecModel.annotation.description. % % The ecModel is saved to the ecModel-specific models/ subfolder. For % saving to other locations, writeYAMLmodel or exportModel are more @@ -91,7 +92,11 @@ function saveEcModel(ecModel,varargin) % model_adapter.toml) and the release that wrote it. ecModel.annotation.gecko_format = '4'; ecModel.annotation.generator = ['GECKO ' geckoVersion()]; - if ~isfield(ecModel.annotation, 'description') || isempty(ecModel.annotation.description) + hasDescription = false; + for f = {'description', 'note'} + hasDescription = hasDescription || (isfield(ecModel.annotation, f{1}) && ~isempty(ecModel.annotation.(f{1}))); + end + if ~hasDescription ecModel.annotation.description = ['Enzyme-constrained model of ' ecModel.id]; end if isfield(ecModel.annotation, 'geckopy_version') diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index a6c227c5..20566514 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -242,7 +242,7 @@ function testsaveECModel_tc0009(testCase) delete(fullfile(adapter.params.path,'models','ecModel.yml')); verifyEqual(testCase, loadedEcModel.annotation.gecko_format, '4') verifyTrue(testCase, startsWith(loadedEcModel.annotation.generator, 'GECKO ')) - loadedEcModel.annotation = rmfield(loadedEcModel.annotation, {'gecko_format', 'generator', 'description'}); + loadedEcModel.annotation = rmfield(loadedEcModel.annotation, {'gecko_format', 'generator', 'note'}); verifyEqual(testCase, ecModel, loadedEcModel) % Test loading of conventional GEM diff --git a/test/unit_tests/parityRound5G4Tests.m b/test/unit_tests/parityRound5G4Tests.m index e6d07c28..61c5aecd 100644 --- a/test/unit_tests/parityRound5G4Tests.m +++ b/test/unit_tests/parityRound5G4Tests.m @@ -50,6 +50,12 @@ function testSaveEcModelWritesDefaultDescriptionOnlyWhenNone_tc0603(testCase) text = fileread(fullfile(project, 'models', 'own.yml')); verifyNotEmpty(testCase, strfind(text, '- description: my own description')) verifyEmpty(testCase, strfind(text, 'Enzyme-constrained model of')) + model.annotation.description = ''; + model.annotation.note = 'loaded description'; + saveEcModel(model, 'note.yml', adapter); + text = fileread(fullfile(project, 'models', 'note.yml')); + verifyNotEmpty(testCase, strfind(text, '- note: loaded description')) + verifyEmpty(testCase, strfind(text, 'Enzyme-constrained model of')) end function model = exchangeToyModel() From 767004709a8c628462e7f159e80267db1e1e8182 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 1 Nov 2026 10:26:09 +0100 Subject: [PATCH 30/31] Test saveEcModel and loadEcModel with an ec reaction without enzymes (tc0700) --- test/VerificationMatrix.tsv | 4 +-- test/unit_tests/geckoCoreFunctionTests.m | 42 ++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 2 deletions(-) diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 253402a0..033ce816 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -80,8 +80,8 @@ model_adapter/adapterParameters.tsv tc0064, tc0065 model_adapter/adapterTemplate.toml tc0068 utilities/* Not explicitly tested - The code has existed for a long time. There are a lot of functions in here, unclear which of them are still useful and functional. io/loadConventionalGEM.m tc0009 -io/loadEcModel.m tc0009, tc0099 -io/saveEcModel.m tc0009, tc0099 +io/loadEcModel.m tc0009, tc0099, tc0700 +io/saveEcModel.m tc0009, tc0099, tc0700 tutorials/tutorial_Human-GEM/HumanGEMAdapter.m tcm0002 - does not test all aspects of the adapter, but some. tutorials/full_ecModel/model_adapter.toml tcm0001 - does not test all aspects of the adapter, but some. tutorials/tutorial_yeast-GEM/code/changeMedia_batch.m Not explicitly tested - The code has existed for a long time diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 4b92695c..2f1ec909 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -252,6 +252,48 @@ function testsaveECModel_tc0009(testCase) verifyEqual(testCase, model, loadedModel) end +function testsaveECModelRxnWithoutEnzymes_tc0700(testCase) + % An ec reaction without enzymes is written without an enzymes entry; + % loadEcModel returns a zero row in rxnEnzMat for it and keeps its id, + % kcat and source. The two forms that earlier versions wrote for it, a + % content-less "enzymes: !!omap" and "enzymes: !!omap []", load the same. + geckoPath = findGECKOroot; + adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m')); + ecModel = makeEcModel(getGeckoTestModel(), false, adapter); + k = find(full(sum(ecModel.ec.rxnEnzMat ~= 0, 2)) > 0, 1); + ecModel.ec.rxnEnzMat(k,:) = 0; + ecModel.ec.source{k} = 'standard'; + file = fullfile(adapter.params.path,'models','ecModel.yml'); + testCase.addTeardown(@() delete(file)); + saveEcModel(ecModel, [], adapter); + lines = strsplit(fileread(file), newline, 'CollapseDelimiters', false); + ecStart = find(strcmp(lines, '- ec-rxns:'), 1); + start = ecStart + find(strcmp(strtrim(lines(ecStart+1:end)), ['- id: ' ecModel.ec.rxns{k}]), 1); + verifyNotEmpty(testCase, start) + % the entry is followed by the next "- !!omap" or by the next section + stop = start + find(~startsWith(lines(start+1:end), ' '), 1); + entry = lines(start:stop-1); + verifyFalse(testCase, any(contains(entry, 'enzymes'))) + verifyTrue(testCase, any(contains(entry, 'source: standard'))) + verifyEqual(testCase, nnz(contains(lines, '!!omap []')), 0) + + forms = {' - enzymes: !!omap', ' - enzymes: !!omap []'}; + for i = 0:numel(forms) + if i > 0 + old = [lines(1:start), forms(i), lines(start+1:end)]; + fid = fopen(file, 'w'); + fprintf(fid, '%s', strjoin(old, newline)); + fclose(fid); + end + loaded = loadEcModel([], adapter); + verifyEqual(testCase, loaded.ec.rxns, ecModel.ec.rxns) + verifyEqual(testCase, loaded.ec.kcat, ecModel.ec.kcat, 'AbsTol', 1e-9) + verifyEqual(testCase, loaded.ec.source, ecModel.ec.source) + verifyEqual(testCase, full(loaded.ec.rxnEnzMat), full(ecModel.ec.rxnEnzMat)) + verifyTrue(testCase, all(loaded.ec.rxnEnzMat(k,:) == 0)) + end +end + function testfuzzyKcatMatchingMedian_tc0060(testCase) % Verifies that kcatAggregation selects which BRENDA aggregate is used. % The fixture records a different median from the max for every triple, From 87017dc23fa66bdc663e189732fd07a53d6556be Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 1 Nov 2026 13:23:49 +0100 Subject: [PATCH 31/31] ecFSEOF orders targets by the rounded score; getECfromGEM warns when the model has no eccodes --- src/flux_analysis/ecFSEOF.m | 6 ++++-- src/get_enzyme_data/getECfromGEM.m | 6 ++++-- test/unit_tests/geckoCoreFunctionTests.m | 9 +++++++++ 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/flux_analysis/ecFSEOF.m b/src/flux_analysis/ecFSEOF.m index df24c4b7..f5801526 100644 --- a/src/flux_analysis/ecFSEOF.m +++ b/src/flux_analysis/ecFSEOF.m @@ -194,8 +194,10 @@ corr_rxns = corr_rxns(targets); target_type = target_type(targets); -% Order from highest to lowest score, the product of slope and correlation -[~,order] = sort(slope_rxns .* corr_rxns,'descend'); +% Order from highest to lowest score, the product of slope and correlation, +% rounded to 9 significant digits so that scores equal up to solver noise keep +% the order of the model +[~,order] = sort(round(slope_rxns .* corr_rxns,9,'significant'),'descend'); target_rxns = target_rxns(order); v_matrix = v_matrix(order,:); rxnGeneM = rxnGeneM(order,:); diff --git a/src/get_enzyme_data/getECfromGEM.m b/src/get_enzyme_data/getECfromGEM.m index dc89922a..c5996fcd 100644 --- a/src/get_enzyme_data/getECfromGEM.m +++ b/src/get_enzyme_data/getECfromGEM.m @@ -3,7 +3,8 @@ % % Uses the model.eccodes to populate the model.ec.eccodes field. EC numbers % must be formatted as four numbers separated by periods, possibly with -% trailing wildcards; incorrectly formatted EC numbers are skipped. +% trailing wildcards; incorrectly formatted EC numbers are skipped. A model without +% a model.eccodes field gives a warning and every EC number is empty. % % Parameters % ---------- @@ -48,7 +49,8 @@ ecRxns = p.ecRxns; if ~isfield(model,'eccodes') - error('The model has no model.eccodes field.') + warning('getECfromGEM:noEccodes','The model has no model.eccodes field, no EC numbers are taken from it.') + model.eccodes = repmat({''},numel(model.rxns),1); end %Need to remove the prefix of GECKO light rxn names in the ec structure diff --git a/test/unit_tests/geckoCoreFunctionTests.m b/test/unit_tests/geckoCoreFunctionTests.m index 98fecda2..6b35794c 100644 --- a/test/unit_tests/geckoCoreFunctionTests.m +++ b/test/unit_tests/geckoCoreFunctionTests.m @@ -184,6 +184,15 @@ function testgetECfromGEM_tc0006(testCase) ecModel2 = getECfromGEM(ecModel,ismember(ecModel.ec.rxns, {'001_R2';'001_R3'})); expEccodes = {'1.1.1.1';'';'1.1.2.1';'';'';''}; verifyEqual(testCase,ecModel2.ec.eccodes,expEccodes) + + %without a model.eccodes field: a warning and empty EC numbers + ecModel = makeEcModel(model, false, adapter); + ecModel = rmfield(ecModel,'eccodes'); + verifyWarning(testCase,@() getECfromGEM(ecModel),'getECfromGEM:noEccodes') + warning('off','getECfromGEM:noEccodes') + ecModel2 = getECfromGEM(ecModel); + warning('on','getECfromGEM:noEccodes') + verifyEqual(testCase,ecModel2.ec.eccodes,repmat({''},6,1)) end %For both full and light