From be3c32bfcca8b4c91ad378c727dd982fd71eb512 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Tue, 3 Nov 2026 12:16:00 +0100 Subject: [PATCH] Measured zero concentrations and isozyme coverage of the proteomics loadProtData reads an abundance of 0 as not measured unless zerosAsMeasured is set. getIsozymeCoverage reports the isozyme groups with a measured and an unmeasured entry, and zeroUnmeasuredIsozymes limits the usage of the unmeasured ones. flexibilizeEnzConcs warns about enzymes with a concentration of 0, which it cannot raise. --- src/io/loadProtData.m | 56 +++- src/limit_proteins/constrainEnzConcs.m | 6 + src/limit_proteins/fillEnzConcs.m | 6 + src/limit_proteins/flexibilizeEnzConcs.m | 30 ++- src/limit_proteins/getIsozymeCoverage.m | 113 ++++++++ src/limit_proteins/levelComplexConcs.m | 22 +- src/limit_proteins/relaxProteomicsGreedy.m | 4 +- src/limit_proteins/zeroUnmeasuredIsozymes.m | 103 +++++++ test/VerificationMatrix.tsv | 10 +- test/unit_tests/isozymeTests.m | 280 ++++++++++++++++++++ test/unit_tests/protDataTests.m | 41 +++ 11 files changed, 647 insertions(+), 24 deletions(-) create mode 100644 src/limit_proteins/getIsozymeCoverage.m create mode 100644 src/limit_proteins/zeroUnmeasuredIsozymes.m create mode 100644 test/unit_tests/isozymeTests.m diff --git a/src/io/loadProtData.m b/src/io/loadProtData.m index 53f10013..36dee3f9 100644 --- a/src/io/loadProtData.m +++ b/src/io/loadProtData.m @@ -24,7 +24,8 @@ % that was previously made by loadProtData. % filterData : logical % whether abundances should be filtered. If false, minVal, maxRSD, -% maxMissing, cutLowest and addStdevs are not considered (default true). +% maxMissing, cutLowest and addStdevs are not considered, and +% zerosAsMeasured applies in both cases. Default true. % modelAdapter : ModelAdapter % a loaded model adapter (default: the current default model adapter). % minVal : double @@ -51,6 +52,18 @@ % how many standard deviations should be added to the mean value of % each protein measurement across replicates, broadening the confidence % interval (default 1). +% zerosAsMeasured : logical +% whether an abundance of 0 is a measurement (default false). If false, +% a 0 is read as not measured (NaN) before any filter, also in a protDataFile +% that is a protData structure: the mean of a protein over its replicates +% leaves out the zeros, and a protein with only zeros is removed. If true, +% a 0 counts as a measured replicate and is kept in the mean; a protein +% whose replicates are all 0 is returned as 0, is not removed by minVal or +% cutLowest, and limits the usage of its enzyme to 0 in constrainEnzConcs. +% Proteomics reports a protein below the detection limit as 0 or as +% missing. With true, a 0 is read as a protein that is absent, for +% example in a deletion strain. maxRSD and maxMissing still apply to a +% protein with some zero replicates. % % Returns % ------- @@ -64,11 +77,16 @@ % - uniprotIDs : cell array with Uniprot IDs matching protData.abundances. % - abundances : matrix of proteomics data, where each column contains mean abundances per condition. % +% A protein that is not measured in a condition is NaN in abundances and +% leaves the usage of its enzyme drawing from the protein pool (see +% fillEnzConcs and constrainEnzConcs). +% % Examples % -------- % % optional arguments may be given positionally or as name-value pairs: -% protData = loadProtData(replPerCond, protDataFile, filterData, modelAdapter, minVal, maxRSD, maxMissing, cutLowest, addStdevs); +% protData = loadProtData(replPerCond, protDataFile, filterData, modelAdapter, minVal, maxRSD, maxMissing, cutLowest, addStdevs, zerosAsMeasured); % protData = loadProtData(replPerCond, 'filterData', false); +% protData = loadProtData(replPerCond, 'zerosAsMeasured', true); % % See also % -------- @@ -82,7 +100,8 @@ 'maxRSD', 1; ... 'maxMissing', 2/3; ... 'cutLowest', 5; ... - 'addStdevs', 1}); + 'addStdevs', 1; ... + 'zerosAsMeasured', false}); protDataFile = p.protDataFile; filterData = p.filterData; modelAdapter = p.modelAdapter; @@ -91,6 +110,10 @@ maxMissing = p.maxMissing; cutLowest = p.cutLowest; addStdevs = p.addStdevs; +zerosAsMeasured = p.zerosAsMeasured; +if isempty(zerosAsMeasured) + zerosAsMeasured = false; +end if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); @@ -129,6 +152,11 @@ abundances = protDataFile.abundances; end +% An abundance of 0 is not measured unless zerosAsMeasured is true +if ~zerosAsMeasured + abundances(abundances == 0) = NaN; +end + %Remove entriew without ID remData = cellfun(@isempty,uniprotIDs); uniprotIDs(remData,:) = []; @@ -149,7 +177,12 @@ else maxMisRepl = maxMissing; end - remData = sum(condAbund>0,2)=0; + else + present = condAbund>0; + end + remData = sum(present,2) 10 + listed = sprintf('%s and %d more', listed, numel(measuredZero) - 10); + end + warning('flexibilizeEnzConcs:measuredZero', ... + ['The enzymes %s have a concentration of 0 in model.ec.concs. Their usage ' ... + 'upper bound stays 0, because it is raised as a factor of the concentration.'], ... + listed) +end % Get enzymes names with a concentration value proteins = model.ec.enzymes(protConcs); @@ -374,7 +394,11 @@ [~, unit] = max(scores); members = units{unit}; - frequence(members(1)) = frequence(members(1)) + 1; + first = members(find(model.ec.concs(protConcs(members)) > 0, 1)); + if isempty(first) + first = members(1); + end + frequence(first) = frequence(first) + 1; unitFrequence(unit) = unitFrequence(unit) + 1; if unitFrequence(unit) > iterPerEnzyme if numel(members) > 1 @@ -404,6 +428,8 @@ end flexIdx = find(frequence > 0 | touched); +% An enzyme with a concentration of 0 keeps its bound and is not reported +flexIdx(model.ec.concs(protConcs(flexIdx)) == 0) = []; protFlex = proteins(flexIdx); oldConcs = model.ec.concs(protConcs(flexIdx)); freq = frequence(flexIdx); diff --git a/src/limit_proteins/getIsozymeCoverage.m b/src/limit_proteins/getIsozymeCoverage.m new file mode 100644 index 00000000..f494b79b --- /dev/null +++ b/src/limit_proteins/getIsozymeCoverage.m @@ -0,0 +1,113 @@ +function [coverage, zeroable] = getIsozymeCoverage(model) +% getIsozymeCoverage Report isozyme groups with partial proteomics coverage. +% +% The reactions in model.ec.rxns that share an identifier once the suffix +% _EXP_ is removed are the isozymes of one reaction of the conventional +% GEM. An entry is measured when at least one of its enzymes has a +% concentration above 0 in model.ec.concs, and unmeasured when none of its +% enzymes has a concentration (all NaN). A group has partial coverage when +% it has both a measured and an unmeasured entry. An unmeasured entry next to +% a measured one can be a protein below the detection limit as well as a +% protein that is absent, so the report makes no statement on which. The kcat +% of an entry is not considered. +% +% Parameters +% ---------- +% model : struct +% an ecModel in full GECKO 3 format (with ecModel.ec structure), not +% GECKO light, with concentrations in model.ec.concs (see fillEnzConcs). +% +% Returns +% ------- +% coverage : table +% one row per group with partial coverage, sorted by reaction, with the +% variables: +% +% - rxn : the reaction identifier without the _EXP_ suffix. +% - nEntries : number of entries of the group in model.ec.rxns. +% - nMeasured : number of measured entries. +% - nUnmeasured : number of unmeasured entries. +% - measuredEnzymes : cell array with the enzymes of the measured entries +% that have a concentration above 0. +% - unmeasuredEnzymes : cell array with the enzymes of the unmeasured +% entries. +% - sharedEnzymes : cell array with the unmeasured enzymes that also +% catalyse an entry outside the unmeasured entries of groups with +% partial coverage. zeroUnmeasuredIsozymes leaves them unchanged. +% zeroable : cell array +% the enzymes that zeroUnmeasuredIsozymes acts on: the enzymes in +% unmeasuredEnzymes that are not in sharedEnzymes of any group. +% +% Examples +% -------- +% coverage = getIsozymeCoverage(model); +% [coverage, zeroable] = getIsozymeCoverage(model); +% +% See also +% -------- +% zeroUnmeasuredIsozymes, getKcatAcrossIsozymes, fillEnzConcs + +if isfield(model, 'ec') && isfield(model.ec, 'geckoLight') && model.ec.geckoLight + error('getIsozymeCoverage:geckoLight', ... + 'This function does not work on GECKO light models.') +end +if ~isfield(model, 'ec') || ~isfield(model.ec, 'rxnEnzMat') || isempty(model.ec.rxnEnzMat) + error('getIsozymeCoverage:noEnzymeMatrix', ... + 'The model has no ec.rxnEnzMat. Run makeEcModel first.') +end + +enzymes = model.ec.enzymes(:); +nEnz = numel(enzymes); +concs = nan(nEnz, 1); +if isfield(model.ec, 'concs') + nConcs = min(numel(model.ec.concs), nEnz); + concs(1:nConcs) = model.ec.concs(1:nConcs); +end +matrix = logical(model.ec.rxnEnzMat); + +entryMeasured = (matrix * double(concs > 0)) > 0; +entryUnmeasured = (matrix * double(~isnan(concs))) == 0; +baseRxn = regexprep(model.ec.rxns, '_EXP_\d+', ''); +[groups, ~, groupIdx] = unique(baseRxn); + +partial = false(numel(groups), 1); +for g = 1:numel(groups) + inGroup = groupIdx == g; + partial(g) = any(entryMeasured(inGroup)) && any(entryUnmeasured(inGroup)); +end + +% An unmeasured entry of a group with partial coverage +inPartial = entryUnmeasured & partial(groupIdx); +% An enzyme is zeroable when it is unmeasured and every entry it catalyses +% is such an entry +catalyses = matrix'; +nEntries = full(sum(catalyses, 2)); +nInPartial = full(catalyses * double(inPartial)); +zeroableMask = isnan(concs) & nEntries > 0 & nInPartial == nEntries; +zeroable = enzymes(zeroableMask); + +rows = find(partial); +n = numel(rows); +rxn = groups(rows); +nEntriesGroup = zeros(n, 1); +nMeasured = zeros(n, 1); +nUnmeasured = zeros(n, 1); +measuredEnzymes = cell(n, 1); +unmeasuredEnzymes = cell(n, 1); +sharedEnzymes = cell(n, 1); +for k = 1:n + inGroup = groupIdx == rows(k); + nEntriesGroup(k) = sum(inGroup); + nMeasured(k) = sum(inGroup & entryMeasured); + nUnmeasured(k) = sum(inGroup & entryUnmeasured); + measuredMask = any(matrix(inGroup & entryMeasured, :), 1)' & concs > 0; + unmeasuredMask = any(matrix(inGroup & entryUnmeasured, :), 1)' & isnan(concs); + measuredEnzymes{k} = enzymes(measuredMask); + unmeasuredEnzymes{k} = enzymes(unmeasuredMask); + sharedEnzymes{k} = enzymes(unmeasuredMask & ~zeroableMask); +end +coverage = table(rxn, nEntriesGroup, nMeasured, nUnmeasured, measuredEnzymes, ... + unmeasuredEnzymes, sharedEnzymes, ... + 'VariableNames', {'rxn','nEntries','nMeasured','nUnmeasured', ... + 'measuredEnzymes','unmeasuredEnzymes','sharedEnzymes'}); +end diff --git a/src/limit_proteins/levelComplexConcs.m b/src/limit_proteins/levelComplexConcs.m index 9b927996..cea82d5d 100644 --- a/src/limit_proteins/levelComplexConcs.m +++ b/src/limit_proteins/levelComplexConcs.m @@ -12,11 +12,13 @@ % Run after fillEnzConcs and before constrainEnzConcs. The concentrations % (mg/gDCW) in model.ec.concs are converted to molar amounts with model.ec.mw % for the levelling and back afterwards. Per complex, the statistic over the -% measured subunits that belong to the complex of the molar amount divided -% by the copies of the subunit is written back, times the copies and the -% molecular weight, to every subunit that belongs to the complex. Subunits +% subunits with a concentration above 0 that belong to the complex of the +% molar amount divided by the copies of the subunit is written back, times +% the copies and the molecular weight, to every subunit that belongs to the +% complex. Subunits % that are shared between complexes, or that also catalyse a reaction on -% their own, are left as measured. +% their own, are left as measured. A subunit with a concentration of 0 is +% measured as absent and stays 0. % % Parameters % ---------- @@ -27,8 +29,8 @@ % Name-Value Arguments % -------------------- % statistic : char -% the summary over the measured subunits, 'median' or 'mean' (default -% 'median'). +% the summary over the subunits with a concentration above 0, 'median' or +% 'mean' (default 'median'). % fillMissing : logical % whether subunits without a measurement, that belong to a complex, are % given the concentration of the complex (default true). @@ -39,7 +41,8 @@ % Returns % ------- % model : struct -% ecModel where model.ec.concs has one abundance per complex. +% ecModel where model.ec.concs has one abundance per complex, except +% for subunits at 0, shared subunits and skipped complexes. % result : struct % what was changed, with the fields: % @@ -51,7 +54,7 @@ % measured because they belong to more than one complex or also % catalyse a reaction on their own. % - skipped : cell array with one cell array of subunits per complex -% that has no measured subunit that belongs to it. +% that has no subunit with a concentration above 0 that belongs to it. % % Examples % -------- @@ -96,6 +99,7 @@ values = molar; isKnown = isfinite(molar); isLevel = isKnown & molar > 0; +isZero = isKnown & molar == 0; result.adjusted = cell(0, 1); result.added = cell(0, 1); result.shared = sort(enzymes(ismember(enzymes, families.shared) & isKnown)); @@ -120,7 +124,7 @@ result.adjusted{end+1, 1} = enzymes{e}; %#ok end values(e) = level; - elseif fillMissing + elseif fillMissing && ~isZero(e) values(e) = level; isKnown(e) = true; result.added{end+1, 1} = enzymes{e}; %#ok diff --git a/src/limit_proteins/relaxProteomicsGreedy.m b/src/limit_proteins/relaxProteomicsGreedy.m index 60171c3c..3d9ad4c9 100644 --- a/src/limit_proteins/relaxProteomicsGreedy.m +++ b/src/limit_proteins/relaxProteomicsGreedy.m @@ -5,7 +5,9 @@ % While growth is below minimalGrowth, finds the still-proteomics-constrained % enzyme (usage_prot_ upper bound below defaultUpperBound) with the % largest absolute shadow price on its prot_ mass-balance constraint, -% fully relaxes it back to defaultUpperBound, re-solves, and repeats. +% fully relaxes it back to defaultUpperBound, re-solves, and repeats. An +% enzyme with a concentration of 0 in model.ec.concs is a candidate like any +% other and is relaxed in the same way. % % Parameters % ---------- diff --git a/src/limit_proteins/zeroUnmeasuredIsozymes.m b/src/limit_proteins/zeroUnmeasuredIsozymes.m new file mode 100644 index 00000000..ef1a9522 --- /dev/null +++ b/src/limit_proteins/zeroUnmeasuredIsozymes.m @@ -0,0 +1,103 @@ +function [model, changed] = zeroUnmeasuredIsozymes(model, varargin) +% zeroUnmeasuredIsozymes Limit unmeasured isozymes that have a measured sibling. +% +% For a reaction with a measured and an unmeasured isozyme (see +% getIsozymeCoverage), the usage upper bound of the enzymes of the unmeasured +% isozyme is lowered, instead of leaving them to draw from the protein pool +% without a limit. The function lowers the usage upper bound in the model +% and does not change model.ec.concs. An enzyme that also catalyses a reaction +% without a measured sibling is left unchanged, because limiting it would +% limit that reaction as well. +% +% Not detecting a protein is not the same as the protein being absent: a +% protein below the detection limit of the measurement is missing from the +% data as well. Method 'zero' treats the unmeasured isozyme as absent and +% sets its usage upper bound to 0, which removes every flux through the +% enzyme and can lower the growth rate or make the model infeasible. Method +% 'detectionLimit' treats it as present at or below the detection limit and +% sets the upper bound to the detection limit. +% +% An enzyme lowered by this function has no concentration in model.ec.concs +% and is a candidate in relaxProteomicsGreedy like any other enzyme with an +% upper bound below the default. +% +% Run after constrainEnzConcs, which resets the usage upper bounds of all +% enzymes. The function is not called by other functions of GECKO. +% +% Parameters +% ---------- +% model : struct +% an ecModel in full GECKO 3 format (with ecModel.ec structure), not +% GECKO light, with concentrations in model.ec.concs, constrained by +% constrainEnzConcs. +% +% Name-Value Arguments +% -------------------- +% method : char +% 'zero' or 'detectionLimit' (default 'zero'). +% limit : double +% the detection limit in mg/gDCW for method 'detectionLimit' (default +% the smallest concentration above 0 in model.ec.concs). +% +% Returns +% ------- +% model : struct +% an ecModel with lowered usage upper bounds. An upper bound that is +% already below the new value is not raised. +% changed : table +% one row per enzyme that the function acted on, with the variables +% uniprot and newUpperBound. +% +% Examples +% -------- +% % optional arguments may be given positionally or as name-value pairs: +% model = zeroUnmeasuredIsozymes(model); +% model = zeroUnmeasuredIsozymes(model, 'method', 'detectionLimit'); +% +% See also +% -------- +% getIsozymeCoverage, constrainEnzConcs, fillEnzConcs + +p = parseGECKOargs(varargin, { ... + 'method', 'zero'; ... + 'limit', []}); +method = p.method; +limit = p.limit; +if isempty(method) + method = 'zero'; +end +if ~ismember(method, {'zero','detectionLimit'}) + error('zeroUnmeasuredIsozymes:method', ... + 'method must be ''zero'' or ''detectionLimit''.') +end + +[~, zeroable] = getIsozymeCoverage(model); + +if strcmp(method, 'zero') + newUb = zeros(numel(zeroable), 1); +else + if isempty(limit) + positive = model.ec.concs(model.ec.concs > 0); + if isempty(positive) + limit = NaN; + else + limit = min(positive); + end + end + if ~isscalar(limit) || ~isnumeric(limit) || isnan(limit) || limit < 0 + error('zeroUnmeasuredIsozymes:limit', ... + 'The detection limit must be a nonnegative number.') + end + newUb = repmat(double(limit), numel(zeroable), 1); +end + +[found, usageIdx] = ismember(strcat('usage_prot_', zeroable), model.rxns); +if any(~found) + error('zeroUnmeasuredIsozymes:noUsageReaction', ... + 'The model has no usage reaction for %s.', strjoin(zeroable(~found), ', ')) +end +newUb = min(newUb(:), model.ub(usageIdx(:))); +usageIdx = usageIdx(:); +model.ub(usageIdx) = newUb; +changed = table(zeroable(:), newUb, 'VariableNames', {'uniprot','newUpperBound'}); +end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 6457127e..345d6623 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -42,13 +42,13 @@ 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 +io/loadProtData.m tc0101, tc0102, tc0260, tc0298, tc0299, tc0313 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/constrainEnzConcs.m tc0503, tc0296, tc0297 limit_proteins/constrainFluxData.m tc0505 limit_proteins/fillEnzConcs.m tc0504 limit_proteins/relaxProteomicsGreedy.m tc0511 @@ -85,9 +85,11 @@ limit_proteins/flexibilizeProteins.m Not explicitly tested - The code has existe limit_proteins/readProteomics.m tc0013 limit_proteins/flexibilizeProtConcs.m tc0013 limit_proteins/getConcControlCoeffs.m tc0013 -limit_proteins/flexibilizeEnzConcs.m tc0013, tc0273, tc0274, tc0275, tc0276 (perComplex option) +limit_proteins/flexibilizeEnzConcs.m tc0013, tc0273, tc0274, tc0275, tc0276 (perComplex option), tc0296, tc0297, tc0317 (measured zero) limit_proteins/getComplexFamilies.m tc0262, tc0263, tc0264, tc0265, tc0266 -limit_proteins/levelComplexConcs.m tc0267, tc0268, tc0269, tc0270, tc0271, tc0272 +limit_proteins/levelComplexConcs.m tc0267, tc0268, tc0269, tc0270, tc0271, tc0272, tc0315, tc0316 +limit_proteins/getIsozymeCoverage.m tc0287, tc0288, tc0289, tc0290, tc0291, tc0314 +limit_proteins/zeroUnmeasuredIsozymes.m tc0292, tc0293, tc0294, tc0295 limit_proteins/updateProtPool.m model_adapter/ModelAdapter.m Base class for adapters, not explicitly tested. model_adapter/ModelAdapterManager.m tc0008, tc0067 diff --git a/test/unit_tests/isozymeTests.m b/test/unit_tests/isozymeTests.m new file mode 100644 index 00000000..afe52e2b --- /dev/null +++ b/test/unit_tests/isozymeTests.m @@ -0,0 +1,280 @@ +%run this test case with the command +%results = runtests('isozymeTests.m') + +function tests = isozymeTests + tests = functiontests(localfunctions); +end + +function testGetIsozymeCoverageReportsPartialGroups_tc0287(testCase) + % W, X and Y have a measured and an unmeasured entry; Z has one entry only. + % P4 also catalyses Z, so zeroing it would limit Z as well. + model = isozymeModel([3; NaN; 2; NaN; NaN]); + [coverage, zeroable] = getIsozymeCoverage(model); + verifyEqual(testCase, coverage.rxn, {'W';'X';'Y'}) + verifyEqual(testCase, coverage.nEntries, [2;2;2]) + verifyEqual(testCase, coverage.nMeasured, [1;1;1]) + verifyEqual(testCase, coverage.nUnmeasured, [1;1;1]) + verifyEqual(testCase, coverage.measuredEnzymes, {{'P1'};{'P1'};{'P3'}}) + verifyEqual(testCase, coverage.unmeasuredEnzymes, {{'P5'};{'P2'};{'P4'}}) + verifyEqual(testCase, coverage.sharedEnzymes, {cell(0,1);cell(0,1);{'P4'}}) + verifyEqual(testCase, zeroable, {'P2';'P5'}) +end + +function testGetIsozymeCoverageCountsAMeasuredSubunitAsMeasuredEntry_tc0288(testCase) + % The entry V has the subunits P1 (measured) and P2 (not measured): it is a + % measured entry, and P2 is not an unmeasured isozyme. + model = isozymeModel([3; NaN; NaN; NaN; NaN]); + model.ec.rxns = {'V';'V_EXP_1'}; + model.ec.rxnEnzMat = sparse([1 1 0 0 0; 0 0 0 1 1]); + [coverage, zeroable] = getIsozymeCoverage(model); + verifyEqual(testCase, coverage.rxn, {'V'}) + verifyEqual(testCase, coverage.measuredEnzymes, {{'P1'}}) + verifyEqual(testCase, coverage.unmeasuredEnzymes, {{'P4';'P5'}}) + verifyEqual(testCase, zeroable, {'P4';'P5'}) +end + +function testGetIsozymeCoverageIgnoresMeasuredZero_tc0289(testCase) + % A concentration of 0 is neither a measured sibling nor an unmeasured entry + model = isozymeModel([0; NaN; 2; NaN; NaN]); + [coverage, zeroable] = getIsozymeCoverage(model); + verifyEqual(testCase, coverage.rxn, {'Y'}) + verifyEmpty(testCase, zeroable) +end + +function testGetIsozymeCoverageEmptyWithoutMeasurements_tc0290(testCase) + model = isozymeModel(nan(5, 1)); + [coverage, zeroable] = getIsozymeCoverage(model); + verifyEqual(testCase, height(coverage), 0) + verifyEmpty(testCase, zeroable) + model.ec.concs = []; + verifyEqual(testCase, height(getIsozymeCoverage(model)), 0) +end + +function testGetIsozymeCoverageRejectsLightModel_tc0291(testCase) + model = isozymeModel([3; NaN; 2; NaN; NaN]); + model.ec.geckoLight = true; + verifyError(testCase, @() getIsozymeCoverage(model), 'getIsozymeCoverage:geckoLight') +end + +function testZeroUnmeasuredIsozymesSetsUsageBoundToZero_tc0292(testCase) + model = isozymeModel([3; NaN; 2; NaN; NaN]); + concs = model.ec.concs; + [model, changed] = zeroUnmeasuredIsozymes(model); + verifyEqual(testCase, changed.uniprot, {'P2';'P5'}) + verifyEqual(testCase, changed.newUpperBound, [0; 0]) + verifyEqual(testCase, usageUb(model, {'P1','P2','P3','P4','P5'}), [1000 0 1000 1000 0]) + verifyEqual(testCase, model.ec.concs, concs) +end + +function testZeroUnmeasuredIsozymesDetectionLimit_tc0293(testCase) + model = isozymeModel([3; NaN; 2; NaN; NaN]); + % The default limit is the smallest concentration above 0 + modelDefault = zeroUnmeasuredIsozymes(model, 'method', 'detectionLimit'); + verifyEqual(testCase, usageUb(modelDefault, {'P2','P5'}), [2 2]) + modelLimit = zeroUnmeasuredIsozymes(model, 'method', 'detectionLimit', 'limit', 0.5); + verifyEqual(testCase, usageUb(modelLimit, {'P2','P5','P4'}), [0.5 0.5 1000]) + % An upper bound below the limit is not raised + model.ub(strcmp(model.rxns, 'usage_prot_P2')) = 0.1; + modelLow = zeroUnmeasuredIsozymes(model, 'method', 'detectionLimit', 'limit', 0.5); + verifyEqual(testCase, usageUb(modelLow, {'P2','P5'}), [0.1 0.5]) +end + +function testZeroUnmeasuredIsozymesNoChangeWithoutPartialGroup_tc0294(testCase) + model = isozymeModel([0; NaN; 2; NaN; NaN]); + ub = model.ub; + [model, changed] = zeroUnmeasuredIsozymes(model); + verifyEqual(testCase, model.ub, ub) + verifyEqual(testCase, height(changed), 0) +end + +function testZeroUnmeasuredIsozymesRejectsUnknownMethod_tc0295(testCase) + model = isozymeModel([3; NaN; 2; NaN; NaN]); + verifyError(testCase, @() zeroUnmeasuredIsozymes(model, 'method', 'half'), ... + 'zeroUnmeasuredIsozymes:method') + verifyError(testCase, @() zeroUnmeasuredIsozymes(model, 'method', 'detectionLimit', 'limit', -1), ... + 'zeroUnmeasuredIsozymes:limit') + noConcs = isozymeModel(nan(5, 1)); + verifyError(testCase, @() zeroUnmeasuredIsozymes(noConcs, 'method', 'detectionLimit'), ... + 'zeroUnmeasuredIsozymes:limit') +end + +function testFlexibilizeEnzConcsLeavesMeasuredZeroAtZero_tc0296(testCase) + % E1 is measured at 1 and limits growth; E2 is measured at 0 on a parallel + % route. E1 is raised, E2 stays at 0 and the function warns about it. + adapter = getAdapter(); + model = zeroToy([1 0]); + verifyWarning(testCase, @() flexibilizeEnzConcs(model, 'expGrowth', 2, 'modelAdapter', adapter, 'verbose', false), ... + 'flexibilizeEnzConcs:measuredZero') + warning('off', 'flexibilizeEnzConcs:measuredZero') + restore = onCleanup(@() warning('on', 'flexibilizeEnzConcs:measuredZero')); + [~, model, flexEnz] = evalc("flexibilizeEnzConcs(model, 'expGrowth', 2, 'modelAdapter', adapter, 'verbose', false)"); + verifyEqual(testCase, flexEnz.uniprotIDs, {'E1'}) + verifyEqual(testCase, usageUb(model, {'E2'}), 0) + sol = solveLP(model); + verifyEqual(testCase, sol.f, 2, 'RelTol', 5e-3) +end + +function testFlexibilizeEnzConcsCannotRaiseMeasuredZero_tc0297(testCase) + % The only route is measured at 0: no enzyme is limiting, the bound stays 0. + % The test pins this behaviour, which a change that raises a bound of 0 breaks. + adapter = getAdapter(); + model = zeroToy([NaN 0]); + model.ub(strcmp(model.rxns, 'R_AB1')) = 0; + warning('off', 'flexibilizeEnzConcs:measuredZero') + restore = onCleanup(@() warning('on', 'flexibilizeEnzConcs:measuredZero')); + [~, model, flexEnz] = evalc("flexibilizeEnzConcs(model, 'expGrowth', 2, 'modelAdapter', adapter, 'verbose', false)"); + verifyEmpty(testCase, flexEnz.uniprotIDs) + verifyEqual(testCase, usageUb(model, {'E2'}), 0) + verifyEqual(testCase, solveLP(model).f, 0, 'AbsTol', 1e-9) +end + + +function testGetIsozymeCoverageKeepsEnzymeOfMeasuredEntry_tc0314(testCase) + % P2 is not measured, and it is in the measured entry X (with P1) as well + % as in the unmeasured entry X_EXP_1, so it is shared and not zeroable. + model = isozymeModel([3; NaN; 2; NaN; NaN]); + model.ec.rxnEnzMat(1, :) = [1 1 0 0 0]; + [coverage, zeroable] = getIsozymeCoverage(model); + verifyEqual(testCase, coverage.sharedEnzymes{strcmp(coverage.rxn, 'X')}, {'P2'}) + verifyEqual(testCase, zeroable, {'P5'}) +end + +function testLevelComplexConcsKeepsMeasuredZero_tc0315(testCase) + % B is measured at 0 in a complex with A and C: it stays 0 and is not + % filled with the level of the complex + model.ec.enzymes = {'A';'B';'C'}; + model.ec.rxns = {'R1'}; + model.ec.genes = {'g_A';'g_B';'g_C'}; + model.ec.mw = [100; 100; 100]; + model.ec.concs = [2; 0; 4]; + model.ec.rxnEnzMat = sparse([1 1 1]); + [model, result] = levelComplexConcs(model); + verifyEqual(testCase, model.ec.concs, [3; 0; 3], 'AbsTol', 1e-12) + verifyEmpty(testCase, result.added) +end + + +function testLevelComplexConcsSkipsComplexWithOnlyZeroSubunits_tc0316(testCase) + model.ec.enzymes = {'A';'B';'C'}; + model.ec.rxns = {'R1'}; + model.ec.genes = {'g_A';'g_B';'g_C'}; + model.ec.mw = [100; 100; 100]; + model.ec.concs = [0; NaN; NaN]; + model.ec.rxnEnzMat = sparse([1 1 1]); + [model, result] = levelComplexConcs(model); + verifyEqual(testCase, model.ec.concs, [0; NaN; NaN]) + verifyEqual(testCase, numel(result.skipped), 1) + verifyEmpty(testCase, result.added) +end + +function testFlexibilizePerComplexKeepsMeasuredZeroOutOfTheReport_tc0317(testCase) + % X1 has a concentration of 0 but an open usage bound, X2 limits the + % complex. X2 is raised and carries the count; X1 is not reported. + adapter = getAdapter(); + model = complexZeroToy([0 1], [1000 1]); + warning('off', 'flexibilizeEnzConcs:measuredZero') + restore = onCleanup(@() warning('on', 'flexibilizeEnzConcs:measuredZero')); + [~, model, flexEnz] = evalc("flexibilizeEnzConcs(model, 'expGrowth', 2, 'modelAdapter', adapter, 'verbose', false, 'perComplex', true)"); + verifyEqual(testCase, flexEnz.uniprotIDs, {'X2'}) + verifyEqual(testCase, flexEnz.frequence, 1) + verifyFalse(testCase, any(isnan(flexEnz.ratioIncr))) + verifyEqual(testCase, usageUb(model, {'X1'}), 1000) +end + + +function adapter = getAdapter() + adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function ub = usageUb(model, enzymes) + [~, idx] = ismember(strcat('usage_prot_', enzymes), model.rxns); + ub = model.ub(idx)'; +end + +function model = isozymeModel(concs) + % Five enzymes with usage reactions (upper bound 1000) and the entries + % X, X_EXP_1: P1, P2; Y, Y_EXP_1: P3, P4; Z: P4; W, W_EXP_1: P5, P1 + enzymes = {'P1';'P2';'P3';'P4';'P5'}; + model.rxns = strcat('usage_prot_', enzymes); + model.ub = 1000 * ones(5, 1); + model.ec.geckoLight = false; + model.ec.enzymes = enzymes; + model.ec.concs = concs(:); + model.ec.rxns = {'X';'X_EXP_1';'Y';'Y_EXP_1';'Z';'W';'W_EXP_1'}; + model.ec.rxnEnzMat = sparse([1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; ... + 0 0 0 1 0; 0 0 0 0 1; 1 0 0 0 0]); +end + +function model = zeroToy(concs) + % A -> B by E1 or by E2 (R_AB1, R_AB2); growth (R4) equals the flux. + mets = {'A_e';'A_c';'B_c';'prot_pool';'prot_E1';'prot_E2'}; + spec = {'EX_A', {'A_e',-1}, -1000, 0; ... + 'TR_A', {'A_e',-1,'A_c',1}, 0, 1000; ... + 'R_AB1', {'A_c',-1,'B_c',1,'prot_E1',-1}, 0, 1000; ... + 'R_AB2', {'A_c',-1,'B_c',1,'prot_E2',-1}, 0, 1000; ... + 'R4', {'B_c',-1}, 0, 1000; ... + 'prot_pool_exchange', {'prot_pool',1}, 0, 1000; ... + 'usage_prot_E1', {'prot_pool',-1,'prot_E1',1}, 0, 1000; ... + 'usage_prot_E2', {'prot_pool',-1,'prot_E2',1}, 0, 1000}; + S = zeros(numel(mets), size(spec, 1)); + for r = 1:size(spec, 1) + entries = spec{r, 2}; + for n = 1:2:numel(entries) + S(strcmp(mets, entries{n}), r) = entries{n+1}; + end + end + model.id = 'toy'; + model.mets = mets; + model.metNames = mets; + model.rxns = spec(:, 1); + model.rxnNames = spec(:, 1); + model.S = sparse(S); + model.lb = cell2mat(spec(:, 3)); + model.ub = cell2mat(spec(:, 4)); + model.c = double(strcmp(model.rxns, 'R4')); + model.b = zeros(numel(mets), 1); + model.rev = double(model.lb < 0); + model.ec.enzymes = {'E1';'E2'}; + model.ec.rxns = {'R_AB1';'R_AB2'}; + model.ec.genes = {'g_E1';'g_E2'}; + model.ec.mw = [100; 100]; + model.ec.concs = concs(:); + model.ec.rxnEnzMat = sparse([1 0; 0 1]); + model = constrainEnzConcs(model); +end + +function model = complexZeroToy(concs, ubs) + % A -> B is catalysed by the complex X1 + X2; growth (R4) equals the flux. + mets = {'A_e';'A_c';'B_c';'prot_pool';'prot_X1';'prot_X2'}; + spec = {'EX_A', {'A_e',-1}, -1000, 0; ... + 'TR_A', {'A_e',-1,'A_c',1}, 0, 1000; ... + 'R_AB', {'A_c',-1,'B_c',1,'prot_X1',-1,'prot_X2',-1}, 0, 1000; ... + 'R4', {'B_c',-1}, 0, 1000; ... + 'prot_pool_exchange', {'prot_pool',1}, 0, 1000; ... + 'usage_prot_X1', {'prot_pool',-1,'prot_X1',1}, 0, ubs(1); ... + 'usage_prot_X2', {'prot_pool',-1,'prot_X2',1}, 0, ubs(2)}; + S = zeros(numel(mets), size(spec, 1)); + for r = 1:size(spec, 1) + entries = spec{r, 2}; + for n = 1:2:numel(entries) + S(strcmp(mets, entries{n}), r) = entries{n+1}; + end + end + model.id = 'toy'; + model.mets = mets; + model.metNames = mets; + model.rxns = spec(:, 1); + model.rxnNames = spec(:, 1); + model.S = sparse(S); + model.lb = cell2mat(spec(:, 3)); + model.ub = cell2mat(spec(:, 4)); + model.c = double(strcmp(model.rxns, 'R4')); + model.b = zeros(numel(mets), 1); + model.rev = double(model.lb < 0); + model.ec.enzymes = {'X1';'X2'}; + model.ec.rxns = {'R_AB'}; + model.ec.genes = {'g_X1';'g_X2'}; + model.ec.mw = [100; 100]; + model.ec.concs = concs(:); + model.ec.rxnEnzMat = sparse([1 1]); +end diff --git a/test/unit_tests/protDataTests.m b/test/unit_tests/protDataTests.m index 2254c766..1e3fabcc 100644 --- a/test/unit_tests/protDataTests.m +++ b/test/unit_tests/protDataTests.m @@ -85,3 +85,44 @@ function testApplyComplexDataStoichiometryByProteinId_tc0103(testCase) end fclose(fid); end + +function testLoadProtDataReadsZeroAsNotMeasuredByDefault_tc0298(testCase) + % Z1 is 0 in every replicate; Z2 is 0 in one replicate. Both are read as + % not measured: Z1 is dropped, and the mean of Z2 is over the two values. + adapter = getTestAdapter(); + file = writeProteomics({'Z1','Z2','Z3'}, [0 0 0; 0 4 6; 5 5 5]); + for filterData = [true false] + protData = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, ... + 'filterData', filterData, 'addStdevs', 0, 'cutLowest', 0); + verifyEqual(testCase, protData.uniprotIDs, {'Z2';'Z3'}) + verifyEqual(testCase, protData.abundances, [5;5], 'AbsTol', 1e-12) + end +end + +function testLoadProtDataKeepsMeasuredZeroWithOption_tc0299(testCase) + % With zerosAsMeasured Z1 stays 0 and the 0 of Z2 enters its mean + adapter = getTestAdapter(); + file = writeProteomics({'Z1','Z2','Z3'}, [0 0 0; 0 4 6; 5 5 5]); + for filterData = [true false] + protData = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, ... + 'filterData', filterData, 'addStdevs', 0, 'cutLowest', 0, 'zerosAsMeasured', true); + verifyEqual(testCase, protData.uniprotIDs, {'Z1';'Z2';'Z3'}) + verifyEqual(testCase, protData.abundances, [0; 10/3; 5], 'AbsTol', 1e-12) + end +end + +function testLoadProtDataMeasuredZeroIsExemptFromMinValAndCutLowest_tc0313(testCase) + % P1 is 0; P2 to P10 are 2 to 10. minVal 2.5 removes P2, then cutLowest + % 20 percent removes one of the eight positive values (P3) and not P1. + adapter = getTestAdapter(); + ids = strcat('P', string(1:10)); + values = repmat([0; (2:10)'], 1, 3); + file = writeProteomics(cellstr(ids), values); + protData = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, ... + 'addStdevs', 0, 'minVal', 2.5, 'cutLowest', 20, 'zerosAsMeasured', true); + verifyEqual(testCase, protData.uniprotIDs, cellstr(strcat('P', string([1 4:10])))') + verifyEqual(testCase, protData.abundances, [0; (4:10)'], 'AbsTol', 1e-12) + protData = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, ... + 'addStdevs', 0, 'minVal', 2.5, 'cutLowest', 20); + verifyEqual(testCase, protData.uniprotIDs, cellstr(strcat('P', string(4:10)))') +end