From 80138aadc6b33ac31db3f0df4f8dd87eea2ae1c7 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 17 Jul 2026 14:24:44 +0200 Subject: [PATCH] test: port invariants asserted by the raven-toolbox suite MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Python suite asserts behaviour on functions RAVEN also has, where the MATLAB tests only assert a type or a count. Ported the invariants that carry over; skipped everything resting on cobra object semantics, pandas or an API RAVEN does not have. manipulation: - convertToIrrev: a reversible reaction with bounds (-500,1000) keeps (0,1000) and its stoichiometry while the _REV copy gets (0,500) and the negated stoichiometry, and inherits the grRule. The existing test asserted only that every rev flag ended up 0, so the bound arithmetic — the only part that can be wrong — was untested. - findDuplicateRxns: a -> b and b -> a group by default and stay separate with ignoreDirection=false. The function's only parameter had no test. - changeGrRules: replace=false ORs onto the existing rule and adds the gene. Only replace=true was covered. - copyToComps: deleteOriginal makes it a move rather than a copy. - mergeModels: metParam decides whether the same metabolite under two ids unifies. The central design decision of the function had no test. conditions: applyCondition had one test, for bounds. Added the prelude reset, metabolite removal, the charge-balance recompute (the most numerically subtle line in the file) and the biomass delta adding to rather than replacing a coefficient. --- testing/function_tests/tConditions.m | 76 ++++++++++++++++ testing/function_tests/tManipulation.m | 118 +++++++++++++++++++++++++ 2 files changed, 194 insertions(+) diff --git a/testing/function_tests/tConditions.m b/testing/function_tests/tConditions.m index 5621f9a0..5d7ce574 100644 --- a/testing/function_tests/tConditions.m +++ b/testing/function_tests/tConditions.m @@ -12,5 +12,81 @@ function applyConditionSetsBounds(testCase) testCase.verifyEqual(m2.lb(idx), 0); end + function applyConditionResetExchangesZeroesUptake(testCase) + % The prelude reopens exchanges to (0, 1000), i.e. no uptake. + m = testCase.model; + [~, exchangeRxns] = getExchangeRxns(m); + m.lb(exchangeRxns) = -37; + m.ub(exchangeRxns) = 42; + condition.prelude.reset_exchanges = 'both'; + evalc('m2 = applyCondition(m, condition);'); + testCase.verifyEqual(unique(m2.lb(exchangeRxns)), 0); + testCase.verifyEqual(unique(m2.ub(exchangeRxns)), 1000); + end + + function applyConditionRemoveMetZeroesCoefficient(testCase) + % remove_mets drops a metabolite from the cofactor pseudoreaction. + m = tConditions.cofactorModel(); + condition.cofactor_pseudoreaction = struct('rxn_id', 'R_cofactor', ... + 'remove_mets', {{struct('met', 'nadh_c')}}); + evalc('m2 = applyCondition(m, condition);'); + testCase.verifyEqual(full(m2.S(strcmp(m2.mets,'nadh_c'), 1)), 0); + end + + function applyConditionRecomputesChargeBalance(testCase) + % charge_balance_met takes whatever coefficient makes the + % pseudoreaction charge neutral, recomputed after the removals. + m = tConditions.cofactorModel(); + condition.cofactor_pseudoreaction = struct('rxn_id', 'R_cofactor', ... + 'remove_mets', {{struct('met', 'nadh_c')}}, ... + 'charge_balance_met', 'h_c'); + evalc('m2 = applyCondition(m, condition);'); + % Left after removing NADH: atp_c (-1, charge -4) and amp_c + % (+1, charge -2), a residual of +2, so H+ must enter at -2. + testCase.verifyEqual(full(m2.S(strcmp(m2.mets,'h_c'), 1)), -2); + % which is to say the reaction ends up charge neutral + idx = find(m2.S(:,1)); + testCase.verifyEqual(sum(full(m2.S(idx,1)) .* m2.metCharges(idx)), 0); + end + + function applyConditionBiomassDeltaAddsToCoefficient(testCase) + % A stoichiometry delta adds to the existing coefficient rather + % than replacing it. + m = tConditions.cofactorModel(); + m.rxns{2} = 'R_biomass'; m.rxnNames{2} = 'R_biomass'; + before = full(m.S(strcmp(m.mets,'atp_c'), 2)); + condition.biomass_stoichiometry_delta = struct('rxn_id', 'R_biomass', ... + 'add', {{struct('met', 'atp_c', 'coef', -0.5)}}); + evalc('m2 = applyCondition(m, condition);'); + testCase.verifyEqual(full(m2.S(strcmp(m2.mets,'atp_c'), 2)), before - 0.5); + end + + end + + methods (Static, Access = private) + + function m = cofactorModel() + % R_cofactor consumes ATP and NADH; R_other consumes ATP. + m = struct(); + m.id = 'cofac'; + m.rxns = {'R_cofactor';'R_other'}; + m.rxnNames = {'R_cofactor';'R_other'}; + m.mets = {'atp_c';'nadh_c';'amp_c';'h_c'}; + m.metNames = {'ATP';'NADH';'AMP';'H+'}; + m.metComps = [1;1;1;1]; + m.metCharges= [-4;-2;-2;1]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + % R_cofactor R_other + m.S = sparse([ -1 -1; % atp_c + -1 0; % nadh_c + 1 0; % amp_c + 0 0]); % h_c + m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0]; + m.b = zeros(4,1); + m.genes = {}; m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0); + m.metFormulas = {'C10';'C21';'C10';'H'}; + end + end end diff --git a/testing/function_tests/tManipulation.m b/testing/function_tests/tManipulation.m index cebf6a1e..b4361edd 100644 --- a/testing/function_tests/tManipulation.m +++ b/testing/function_tests/tManipulation.m @@ -113,6 +113,86 @@ function convertToIrrevAllIrreversible(testCase) testCase.verifyGreaterThanOrEqual(numel(m2.rxns), numel(testCase.model.rxns)); end + function convertToIrrevSplitsBoundsAndStoichiometry(testCase) + % A reversible reaction with bounds (-500,1000) keeps (0,1000) and + % its stoichiometry, while the _REV copy gets (0,500) and the + % negated stoichiometry. + m = tManipulation.twoMetModel(); + m.rev = 1; m.lb = -500; m.ub = 1000; + m2 = convertToIrrev(m); + + fwd = strcmp(m2.rxns, 'R1'); + rev = strcmp(m2.rxns, 'R1_REV'); + testCase.verifyTrue(any(rev)); + testCase.verifyEqual(full(m2.lb(fwd)), 0); + testCase.verifyEqual(full(m2.ub(fwd)), 1000); + testCase.verifyEqual(full(m2.lb(rev)), 0); + testCase.verifyEqual(full(m2.ub(rev)), 500); + testCase.verifyEqual(full(m2.S(:,fwd)), [-1;1]); + testCase.verifyEqual(full(m2.S(:,rev)), [1;-1]); + end + + function convertToIrrevReverseCopyInheritsGrRule(testCase) + % The _REV copy is catalysed by the same genes as the forward one. + m = tManipulation.twoMetModel(); + m.rev = 1; m.lb = -500; m.ub = 1000; + m.genes = {'g1'}; m.grRules = {'g1'}; m.rxnGeneMat = sparse(1,1,1); + m2 = convertToIrrev(m); + testCase.verifyEqual(m2.grRules{strcmp(m2.rxns,'R1_REV')}, 'g1'); + end + + function findDuplicateRxnsIgnoreDirection(testCase) + % a -> b and b -> a are the same reaction run backwards, so they + % group by default and stay separate when direction matters. + m = tManipulation.twoMetModel(); + m.rxns = {'R1';'R2'}; + m.rxnNames = {'R1';'R2'}; + m.S = sparse([-1 1; 1 -1]); + m.lb = [0;0]; m.ub = [1000;1000]; m.rev = [0;0]; m.c = [0;0]; + m.grRules = {'';''}; m.rxnGeneMat = sparse(2,0); + + pairs = findDuplicateRxns(m); + testCase.verifyEqual(pairs, [1 2]); + + pairs = findDuplicateRxns(m, 'ignoreDirection', false); + testCase.verifyEmpty(pairs); + end + + function changeGrRulesAppendsToExistingRule(testCase) + % replace=false must OR the new rule onto the existing one rather + % than overwrite it, and add the new gene to the model. + m2 = changeGrRules(testCase.model, 'ACKr', 'b9999', false); + rule = m2.grRules{strcmp(m2.rxns, 'ACKr')}; + testCase.verifySubstring(rule, 'b9999'); + testCase.verifySubstring(rule, ' or '); + % the gene the reaction already had must still be in the rule + oldRule = testCase.model.grRules{strcmp(testCase.model.rxns, 'ACKr')}; + oldGene = regexp(oldRule, 'b\d+', 'match', 'once'); + testCase.verifySubstring(rule, oldGene); + testCase.verifyTrue(ismember('b9999', m2.genes)); + end + + function copyToCompsDeleteOriginalIsAMove(testCase) + % deleteOriginal turns the copy into a move: the reaction count is + % unchanged and the original compartment's copy is gone. + evalc('copied = copyToComps(testCase.model, {''p''}, ''rxns'', ''ACKr'');'); + evalc(['moved = copyToComps(testCase.model, {''p''}, ''rxns'', ''ACKr'', ' ... + '''deleteOriginal'', true);']); + testCase.verifyEqual(numel(copied.rxns), numel(testCase.model.rxns) + 1); + testCase.verifyEqual(numel(moved.rxns), numel(testCase.model.rxns)); + end + + function mergeModelsMetParamDecidesUnification(testCase) + % The same metabolite under two ids unifies when matching on names + % and stays distinct when matching on ids. + a = tManipulation.namedMetModel('glc_c', 'A'); + b = tManipulation.namedMetModel('glucose_c', 'B'); + evalc('byName = mergeModels({a; b});'); + evalc('byId = mergeModels({a; b}, ''metParam'', ''mets'');'); + testCase.verifyEqual(nnz(strcmp(byName.metNames, 'Glucose')), 1); + testCase.verifyEqual(nnz(strcmp(byId.metNames, 'Glucose')), 2); + end + function copyToCompsAddsCompartment(testCase) evalc('m2 = copyToComps(testCase.model, {''p''}, ''ACKr'');'); testCase.verifyTrue(ismember('p', m2.comps)); @@ -252,4 +332,42 @@ function standardizeGrRulesReturnsRules(testCase) end end + + methods (Static, Access = private) + + function m = twoMetModel() + % Single reaction a -> b in one compartment. + m = struct(); + m.id = 'toy'; + m.rxns = {'R1'}; + m.rxnNames = {'R1'}; + m.mets = {'a';'b'}; + m.metNames = {'a';'b'}; + m.metComps = [1;1]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.S = sparse([-1;1]); + m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; m.b = zeros(2,1); + m.genes = {}; m.grRules = {''}; m.rxnGeneMat = sparse(1,0); + m.metFormulas = {'C';'C'}; + end + + function m = namedMetModel(glucoseId, modelId) + % Glucose[c] under a caller-chosen id, consumed by one reaction. + m = struct(); + m.id = modelId; + m.rxns = {['R_' modelId]}; + m.rxnNames = m.rxns; + m.mets = {glucoseId; ['product_' modelId]}; + m.metNames = {'Glucose'; ['Product' modelId]}; + m.metComps = [1;1]; + m.comps = {'c'}; + m.compNames = {'cytosol'}; + m.S = sparse([-1;1]); + m.lb = 0; m.ub = 1000; m.rev = 0; m.c = 0; m.b = zeros(2,1); + m.genes = {}; m.grRules = {''}; m.rxnGeneMat = sparse(1,0); + m.metFormulas = {'C6H12O6';'C'}; + end + + end end