Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions src/enzyme_usage/getEnzymeUsage.m
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,25 @@
% usage_prot_<id> reaction in the model is left out, with the warning
% getEnzymeUsage:missingUsageReaction.
% - fluxes : vector of fluxes, copy of input fluxes.
% - measuredProtein : sum of the measured concentrations in model.ec.concs
% (mg/gDCW), over the enzymes with a non-NaN entry.
% - demandedProtein : sum of absUsage over the same enzymes (mg/gDCW).
% - demandRatio : demandedProtein / measuredProtein. NaN when no
% concentration is measured; Inf when the measured concentrations sum to 0
% and the demand is positive. A ratio above 1 means that the solution uses
% more of the measured enzymes than was measured, which happens when the
% upper bounds of the usage reactions are above the measurements (for
% instance after flexibilizeEnzConcs). A ratio at or below 1 does not
% exclude that one enzyme exceeds its measurement, since the sum is over
% all measured enzymes; capUsage gives the ratio per enzyme. Enzymes
% without a measured concentration draw from the protein pool and are not
% counted.
% - measuredProteinActive, demandRatioActive : the same two values restricted
% to the enzymes with an absUsage above 1e-7 mg/gDCW. Measured enzymes
% that are inactive in the solution raise measuredProtein and can hide a
% conflict among the active ones.
% The five fields are computed over all enzymes with a usage reaction,
% whatever the value of zero.
%
% Examples
% --------
Expand All @@ -61,13 +80,17 @@
error('This function does not work on GECKO light models.')
end
usageData.protID = ecModel.ec.enzymes;
concs = nan(numel(ecModel.ec.enzymes),1);
nConcs = min(numel(ecModel.ec.concs),numel(concs));
concs(1:nConcs) = ecModel.ec.concs(1:nConcs);
[~,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);
concs = concs(~noUsageRxn);
rxnIdx = rxnIdx(~noUsageRxn);
end

Expand All @@ -76,6 +99,14 @@
usageData.capUsage = abs(usageData.absUsage./usageData.UB);
usageData.fluxes = fluxes;

measured = ~isnan(concs);
active = usageData.absUsage > 1e-7;
usageData.measuredProtein = sum(concs(measured));
usageData.measuredProteinActive = sum(concs(measured & active));
usageData.demandedProtein = sum(usageData.absUsage(measured));
usageData.demandRatio = usageData.demandedProtein / usageData.measuredProtein;
usageData.demandRatioActive = sum(usageData.absUsage(measured & active)) / usageData.measuredProteinActive;

if ~zero
nonzero = usageData.absUsage>0;
usageData.absUsage = usageData.absUsage(nonzero);
Expand Down
14 changes: 14 additions & 0 deletions src/enzyme_usage/reportEnzymeUsage.m
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,14 @@
% the topAbsUsage enzymes with the highest absolute usage.
% totalUsageFlux : double
% total protein pool used, i.e. the upper bound of prot_pool_exchange.
% demandRatio : double
% demandRatio of usageData: the protein used of the enzymes with a
% measured concentration divided by their measured concentrations.
% Above 1, the solution needs more of those enzymes than was
% measured. NaN when usageData has no demandRatio field.
% demandRatioActive : double
% demandRatioActive of usageData: the same ratio over the enzymes
% that carry flux in the solution.
%
% Examples
% --------
Expand Down Expand Up @@ -197,4 +205,10 @@
end
usageReport.topAbsUsage = struct2table(topUsage);
usageReport.totalUsageFlux = protPool;
usageReport.demandRatio = NaN;
usageReport.demandRatioActive = NaN;
if isfield(usageData, 'demandRatio')
usageReport.demandRatio = usageData.demandRatio;
usageReport.demandRatioActive = usageData.demandRatioActive;
end
end
4 changes: 2 additions & 2 deletions test/VerificationMatrix.tsv
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,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
enzyme_usage/getEnzymeUsage.m tc0286, tc0507, tc0508
enzyme_usage/reportEnzymeUsage.m tc0286, tc0508, tc0509
gather_kcats/fuzzyKcatMatching.m tc0010, tc0223, tc0224, tc0225, tc0226, tc0231, tc0239
gather_kcats/getStandardKcat.m tc0256, tc0257, tc0309, tc0310
gather_kcats/mergeDlkcatAndFuzzyKcats.m tc0011; tc0241
Expand Down
80 changes: 80 additions & 0 deletions test/unit_tests/geckoCoreFunctionTests.m
Original file line number Diff line number Diff line change
Expand Up @@ -1291,6 +1291,86 @@ function testPfbaEnzymesFractionOfNegativeOptimum_tc0285(testCase)
verifyGreaterThanOrEqual(testCase, solution.objectiveValue, -11)
end

function testGetEnzymeUsageDemandRatio_tc0286(testCase)
% Verifies the demand ratio of getEnzymeUsage: the usage of the enzymes
% with a measured concentration divided by their measured
% concentrations, over all of them and over those that carry flux, NaN
% without a measured concentration, and the same whatever 'zero' is.
geckoPath = findGECKOroot;
adapter = ModelAdapterManager.getAdapter(fullfile(geckoPath,'test','unit_tests','ecTestGEM', 'TestGEMAdapter.m'));
model = getGeckoTestModel();
ecModel = makeEcModel(model, false, adapter);
ecModel = getECfromGEM(ecModel);
ecModel.ec.kcat(:) = 100;
ecModel.ec.source(:) = {'manual'};
ecModel = applyKcatConstraints(ecModel);
ecModel = setProtPoolSize(ecModel, [], [], [], adapter);

[~, idx] = ismember(strcat('usage_prot_', ecModel.ec.enzymes), ecModel.rxns);
ecModel.ec.concs = nan(numel(ecModel.ec.enzymes), 1);
ecModel.ec.concs(1:3) = [2; 4; 1];
fluxes = zeros(numel(ecModel.rxns), 1);
fluxes(idx(1)) = 3;
fluxes(idx(3)) = 1.5;
fluxes(idx(4)) = 7;

usageData = getEnzymeUsage(ecModel, fluxes);
verifyEqual(testCase, usageData.measuredProtein, 7, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.demandedProtein, 4.5, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.demandRatio, 4.5/7, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.measuredProteinActive, 3, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.demandRatioActive, 1.5, 'AbsTol', 1e-12)

withoutZero = getEnzymeUsage(ecModel, fluxes, 'zero', false);
verifyEqual(testCase, withoutZero.demandRatio, usageData.demandRatio, 'AbsTol', 1e-12)
verifyEqual(testCase, withoutZero.demandRatioActive, usageData.demandRatioActive, 'AbsTol', 1e-12)

report = reportEnzymeUsage(ecModel, usageData);
verifyEqual(testCase, report.demandRatio, usageData.demandRatio, 'AbsTol', 1e-12)
verifyEqual(testCase, report.demandRatioActive, usageData.demandRatioActive, 'AbsTol', 1e-12)

% A usage below the threshold counts as demand but not as active.
lowFluxes = fluxes;
lowFluxes(idx(2)) = 1e-8;
usageData = getEnzymeUsage(ecModel, lowFluxes);
verifyEqual(testCase, usageData.demandedProtein, 4.5 + 1e-8, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.measuredProteinActive, 3, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.demandRatioActive, 1.5, 'AbsTol', 1e-12)

% An enzyme without a usage reaction is left out with its concentration.
reduced = removeReactions(ecModel, ecModel.rxns(idx(1)));
reducedFluxes = fluxes;
reducedFluxes(idx(1)) = [];
warning('off', 'getEnzymeUsage:missingUsageReaction')
restore = onCleanup(@() warning('on', 'getEnzymeUsage:missingUsageReaction'));
usageData = getEnzymeUsage(reduced, reducedFluxes);
verifyEqual(testCase, usageData.measuredProtein, 5, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.demandedProtein, 1.5, 'AbsTol', 1e-12)
verifyEqual(testCase, usageData.demandRatioActive, 1.5, 'AbsTol', 1e-12)

% A usageData without the demand fields gives NaN in the report.
bare = rmfield(getEnzymeUsage(ecModel, fluxes), ...
{'measuredProtein','measuredProteinActive','demandedProtein','demandRatio','demandRatioActive'});
report = reportEnzymeUsage(ecModel, bare);
verifyTrue(testCase, isnan(report.demandRatio))
verifyTrue(testCase, isnan(report.demandRatioActive))

% A concentration vector shorter than the enzyme list counts as unmeasured.
short = ecModel;
short.ec.concs = [2; 4];
usageData = getEnzymeUsage(short, fluxes);
verifyEqual(testCase, usageData.measuredProtein, 6, 'AbsTol', 1e-12)
short.ec.concs = [];
usageData = getEnzymeUsage(short, fluxes);
verifyTrue(testCase, isnan(usageData.demandRatio))

ecModel.ec.concs(:) = NaN;
usageData = getEnzymeUsage(ecModel, fluxes);
verifyEqual(testCase, usageData.measuredProtein, 0)
verifyTrue(testCase, isnan(usageData.demandRatio))
verifyTrue(testCase, isnan(usageData.demandRatioActive))
end


function ecModel = fixtureForRelaxProteomicsGreedy()
% Builds the shared fixture used by the relaxProteomicsGreedy tests below:
Expand Down
Loading