diff --git a/src/limit_proteins/flexibilizeEnzConcs.m b/src/limit_proteins/flexibilizeEnzConcs.m index 2e3cf062..b94e5d57 100644 --- a/src/limit_proteins/flexibilizeEnzConcs.m +++ b/src/limit_proteins/flexibilizeEnzConcs.m @@ -39,6 +39,17 @@ % a loaded model adapter (default: the current default model adapter). % verbose : logical % whether progress should be reported (default true). +% perComplex : logical +% whether every complex (see getComplexFamilies) is flexibilized as one +% enzyme (default false). Enzymes and complexes are then ranked by the +% growth gained per foldChange of raising all their usage upper bounds +% together (a levelled complex binds on all subunits at once, so no +% single subunit shows a gain), every measured subunit of the picked +% complex is raised by the same factor, iterations are counted per +% complex, and one factor per complex is kept in the final pass, so that +% concentrations levelled by levelComplexConcs stay levelled. An upper +% bound is never lowered below its current value. The frequence of a +% complex is recorded on its first measured subunit. % % Returns % ------- @@ -66,23 +77,28 @@ % % See also % -------- -% constrainEnzConcs, getConcControlCoeffs +% constrainEnzConcs, getConcControlCoeffs, levelComplexConcs, getComplexFamilies p = parseGECKOargs(varargin, { ... 'expGrowth', []; ... 'foldChange', []; ... 'iterPerEnzyme', []; ... 'modelAdapter', []; ... - 'verbose', []}); + 'verbose', []; ... + 'perComplex', []}); expGrowth = p.expGrowth; foldChange = p.foldChange; iterPerEnzyme = p.iterPerEnzyme; modelAdapter = p.modelAdapter; verbose = p.verbose; +perComplex = p.perComplex; if isempty(verbose) verbose = true; end +if isempty(perComplex) + perComplex = false; +end if isempty(modelAdapter) modelAdapter = ModelAdapterManager.getDefault(); if isempty(modelAdapter) @@ -134,6 +150,12 @@ % Store how many time is predicted a enzyme with the highest coeff frequence = zeros(numel(proteins),1); +if perComplex + [model, flexEnz] = flexibilizePerComplex(model, proteins, protConcs, ... + predGrowth, expGrowth, foldChange, iterPerEnzyme, params, verbose); + return +end + flexBreak=false; if any(protConcs) while predGrowth < expGrowth @@ -265,3 +287,249 @@ flexEnz.flexConcs = flexEnz.flexConcs(b); flexEnz.frequence = flexEnz.frequence(b); end + +function [model, flexEnz] = flexibilizePerComplex(model, proteins, protConcs, predGrowth, expGrowth, foldChange, iterPerEnzyme, params, verbose) +% Flexibilize with every complex as one unit, see the perComplex option. +if isempty(protConcs) + error('Protein concentrations have not been defined. Please run readProteomics and constrainEnzConcs') +end +limit = 0.75; +nProt = numel(proteins); +[~, usageIdx] = ismember(strcat('usage_prot_', proteins), model.rxns); +bioRxnIdx = getIndexes(model, params.bioRxn, 'rxns'); +families = getComplexFamilies(model); + +% Positions in proteins of each unit that is relaxed together: the measured +% subunits that belong to a complex, or a single enzyme +units = {}; +assigned = false(nProt, 1); +for i = 1:nProt + if assigned(i) + continue + end + [~, e] = ismember(proteins{i}, model.ec.enzymes); + c = families.owner(e); + members = i; + if c > 0 + cplx = families.complexes{c}; + [~, members] = ismember(cplx.enzymes(cplx.owned), proteins); + members = members(members > 0)'; + end + assigned(members) = true; + units{end+1, 1} = members; %#ok +end +unitFrequence = zeros(numel(units), 1); +frequence = zeros(nProt, 1); +touched = false(nProt, 1); + +% Backstop on the iterations when the number per unit is not limited +if isinf(iterPerEnzyme) + capPerEnzyme = 1000; +else + capPerEnzyme = iterPerEnzyme; +end +maxTotalIter = nProt * capPerEnzyme + nProt; +iterCount = 0; +flexBreak = false; +changedProtPool = false; + +while predGrowth < expGrowth + iterCount = iterCount + 1; + if iterCount > maxTotalIter + printOrange(['WARNING: The global iteration limit has been reached without reaching expGrowth. ' ... + 'The target may be unreachable even with relaxed enzyme bounds.\n']); + flexBreak = true; + break + end + scores = getUnitScores(model, usageIdx, units, foldChange, limit); + + if ~any(scores > 0) + disp('No (more) limiting enzymes have been found. Attempting to increase protein pool exchange...') + protPoolIdx = getIndexes(model,'prot_pool_exchange','rxns'); + oldProtPool = model.ub(protPoolIdx); + tempModel = setParam(model,'ub',protPoolIdx,1000); + tempModel = setParam(tempModel,'ub',bioRxnIdx,expGrowth); + sol = solveLP(tempModel); + if (sol.f-predGrowth)>1e-3 + predGrowth = sol.f; + tempModel = setParam(tempModel,'lb',bioRxnIdx,predGrowth); + tempModel = setParam(tempModel,'obj',protPoolIdx,-1); + sol = solveLP(tempModel); + newProtPool = sol.x(protPoolIdx); + model.ub(protPoolIdx) = newProtPool; + fprintf(['Changing the upper bound of protein pool exchange from %s ',... + 'to %s enabled a\ngrowth rate of %s. It can be helpful to set ', ... + 'the upper bound of protein\npool exchange to 1000 before ',... + 'running flexibilizeEnzConcs.\n'],num2str(oldProtPool), ... + num2str(newProtPool),num2str(predGrowth)); + changedProtPool = true; + else + fprintf(['Protein pool exchange was also not limiting. Inability to reach growth ',... + 'rate is not related to\nenzyme constraints. Maximum growth rate ',... + 'is %s.\n'], num2str(predGrowth)) + end + expGrowth = predGrowth; + break + end + + [~, unit] = max(scores); + members = units{unit}; + frequence(members(1)) = frequence(members(1)) + 1; + unitFrequence(unit) = unitFrequence(unit) + 1; + if unitFrequence(unit) > iterPerEnzyme + if numel(members) > 1 + what = ['complex of ' strjoin(proteins(members)', ', ')]; + else + what = ['Protein ' proteins{members(1)}]; + end + printOrange(['WARNING: Limit has been reached. ' what ' seems to be problematic. Consider changing the kcat.\n']); + flexBreak = true; + break + end + + increase = foldChange * unitFrequence(unit); + for i = members + touched(i) = true; + model.ub(usageIdx(i)) = max(model.ub(usageIdx(i)), model.ec.concs(protConcs(i)) * (1+increase)); + end + sol = solveLP(model); + if sol.stat == 1 + predGrowth = sol.f; + else + predGrowth = 0; + end + if verbose + disp(['Protein ' proteins{members(1)} ' UB adjusted. Grow: ' num2str(predGrowth)]) + end +end + +flexIdx = find(frequence > 0 | touched); +protFlex = proteins(flexIdx); +oldConcs = model.ec.concs(protConcs(flexIdx)); +freq = frequence(flexIdx); +if isempty(flexIdx) + flexConcs = zeros(0, 1); +elseif flexBreak + flexConcs = model.ub(usageIdx(flexIdx)); +else + % Not all flexibilized proteins require flexibilization in the end. + % Test flux distribution with minimum prot_pool usage + modelTemp = setParam(model,'var',params.bioRxn,expGrowth,0.5); + modelTemp = setParam(modelTemp,'obj','prot_pool_exchange',-1); + sol = solveLP(modelTemp,1); + newConcs = abs(sol.x(usageIdx(flexIdx))); + % One factor per complex: every flexibilized subunit takes the largest + % new/old of its complex + seen = false(numel(flexIdx), 1); + for k = 1:numel(flexIdx) + if seen(k) + continue + end + [~, e] = ismember(protFlex{k}, model.ec.enzymes); + c = families.owner(e); + group = k; + if c > 0 + cplx = families.complexes{c}; + [~, group] = ismember(cplx.enzymes(cplx.owned), protFlex); + group = group(group > 0)'; + end + seen(group) = true; + if numel(group) < 2 + continue + end + withOld = group(oldConcs(group) > 0); + if isempty(withOld) + continue + end + factor = max(newConcs(withOld) ./ oldConcs(withOld)); + newConcs(withOld) = oldConcs(withOld) * factor; + end + keepOld = newConcs < oldConcs; + model.ub(usageIdx(flexIdx(keepOld))) = oldConcs(keepOld); + model.ub(usageIdx(flexIdx(~keepOld))) = newConcs(~keepOld); + protFlex(keepOld) = []; + oldConcs(keepOld) = []; + freq(keepOld) = []; + flexConcs = newConcs(~keepOld); +end + +flexEnz.uniprotIDs = protFlex; +flexEnz.oldConcs = oldConcs; +flexEnz.flexConcs = flexConcs; +flexEnz.frequence = freq; +if changedProtPool + flexEnz.uniprotIDs(end+1) = {'prot_pool'}; + flexEnz.oldConcs(end+1) = oldProtPool; + flexEnz.flexConcs(end+1) = newProtPool; + flexEnz.frequence(end+1) = 1; +end +flexEnz.ratioIncr = flexEnz.flexConcs./flexEnz.oldConcs; +[~,b] = sort(flexEnz.ratioIncr,'descend'); +flexEnz.ratioIncr = flexEnz.ratioIncr(b); +flexEnz.uniprotIDs = flexEnz.uniprotIDs(b); +flexEnz.oldConcs = flexEnz.oldConcs(b); +flexEnz.flexConcs = flexEnz.flexConcs(b); +flexEnz.frequence = flexEnz.frequence(b); +end + +function scores = getUnitScores(model, usageIdx, units, foldChange, limit) +% Growth gained per fold of raising the usage upper bounds of each unit +% together. With reduced costs, the sum over the unit's subunits that run +% above limit of their cap of |reduced cost of the usage reaction| times its +% upper bound. Unlike the shadow price of the protein metabolite, the +% reduced cost of the usage reaction leaves out the price of the protein +% pool, which would otherwise be counted once per subunit. When no such +% score is positive (no reduced costs, or the protein pool binds and the +% solver puts the whole price on it), the candidate units are probed one by +% one with all their upper bounds raised by foldChange and the pool lifted. +scores = zeros(numel(units), 1); +sol = solveLP(model); +if sol.stat ~= 1 + return +end +candidates = zeros(0, 1); +active = cell(numel(units), 1); +for k = 1:numel(units) + rx = usageIdx(units{k}); + ub = model.ub(rx); + act = ub > 0; + act(act) = sol.x(rx(act)) ./ ub(act) > limit; + if any(act) + candidates(end+1, 1) = k; %#ok + active{k} = rx(act); + end +end +if isempty(candidates) + return +end +needed = vertcat(active{candidates}); +if ~isempty(sol.rCost) && all(isfinite(sol.rCost(needed))) + for k = candidates' + scores(k) = sum(abs(sol.rCost(active{k})) .* model.ub(active{k})); + end + scores(scores < 1e-12) = 0; + if any(scores > 0) + return + end + scores = zeros(numel(units), 1); +end +poolIdx = getIndexes(model, 'prot_pool_exchange', 'rxns'); +model.ub(poolIdx) = max(model.ub(poolIdx), 1000); +base = solveLP(model); +if base.stat ~= 1 + return +end +for k = candidates' + tempModel = model; + rx = usageIdx(units{k}); + tempModel.ub(rx) = tempModel.ub(rx) * foldChange; + probe = solveLP(tempModel); + if probe.stat == 1 + gain = probe.f - base.f; + if gain > 1e-3 + scores(k) = gain / (foldChange - 1); + end + end +end +end + diff --git a/src/limit_proteins/getComplexFamilies.m b/src/limit_proteins/getComplexFamilies.m new file mode 100644 index 00000000..f7a27c6c --- /dev/null +++ b/src/limit_proteins/getComplexFamilies.m @@ -0,0 +1,198 @@ +function families = getComplexFamilies(model, varargin) +% getComplexFamilies Group the enzymes of multi-subunit reactions into complexes. +% +% A complex is a set of enzymes that catalyses a reaction together, taken +% from model.ec.rxnEnzMat. Two enzyme sets are one complex when they share +% at least minShared subunits and at least a fraction minOverlap of the +% larger set (isozyme variants that swap a subunit, the same complex on +% different reactions). A single shared protein never merges two sets, so +% electron carriers and shared subunits (cytochrome c, thioredoxin, the E3 +% subunit of the 2-oxo-acid dehydrogenases) do not chain unrelated +% complexes together. +% +% A subunit in exactly two complexes belongs to the one whose enzyme sets +% catalyse more reactions (the core subunits COX1-3 of complex IV also +% appear in one assembly reaction). A subunit in three or more complexes, +% in two with equal weight, or that also catalyses a reaction on its own is +% shared: it is left as measured and does not take part in the summary of a +% complex. +% +% When isozyme variants of a complex disagree on the number of copies of a +% subunit, the first variant is used and a warning is given. +% +% Parameters +% ---------- +% model : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure), with +% ec.rxnEnzMat (copies per subunit). +% +% Name-Value Arguments +% -------------------- +% minOverlap : double +% two enzyme sets are one complex when they share at least this +% fraction of the larger set's subunits (default 0.5)... +% minShared : double +% ...and at least this many subunits (default 2). +% +% Returns +% ------- +% families : struct +% the complexes of the model, with the fields: +% +% - complexes : cell array with one struct per complex, with the fields +% enzymes (cell array with every subunit that appears in any of its +% enzyme sets), copies (copies of each of these subunits) and owned +% (logical, true for the subunits that belong to this complex). +% - owner : numeric vector with one entry per model.ec.enzymes, the +% index in complexes of the complex that the subunit belongs to, or 0 +% if it belongs to none. +% - shared : cell array with the subunits of a complex that belong to +% none: carriers and subunits shared between complexes. +% +% Examples +% -------- +% % optional arguments may be given positionally or as name-value pairs: +% families = getComplexFamilies(model); +% families = getComplexFamilies(model, 'minShared', 3); +% +% See also +% -------- +% levelComplexConcs, flexibilizeEnzConcs + +p = parseGECKOargs(varargin, { ... + 'minOverlap', []; ... + 'minShared', []}); +minOverlap = p.minOverlap; +minShared = p.minShared; +if isempty(minOverlap); minOverlap = 0.5; end +if isempty(minShared); minShared = 2; end + +if ~isfield(model, 'ec') || ~isfield(model.ec, 'rxnEnzMat') || isempty(model.ec.rxnEnzMat) + error('getComplexFamilies:noRxnEnzMat', ... + 'The model has no ec.rxnEnzMat. Run makeEcModel first.') +end +enzymes = model.ec.enzymes; +matrix = sparse(model.ec.rxnEnzMat); + +% Distinct enzyme sets of the reactions with more than one enzyme, in order +% of first appearance, with the copies of the first reaction that has the set +setIdx = {}; +setCopies = {}; +setCount = zeros(0, 1); +setLookup = containers.Map('KeyType', 'char', 'ValueType', 'double'); +alone = false(numel(enzymes), 1); +for row = 1:size(matrix, 1) + [~, cols, vals] = find(matrix(row, :)); + if numel(cols) == 1 + alone(cols) = true; + elseif ~isempty(cols) + key = sprintf('%d,', cols); + if isKey(setLookup, key) + k = setLookup(key); + else + k = numel(setIdx) + 1; + setLookup(key) = k; + setIdx{k, 1} = cols(:); + setCopies{k, 1} = full(vals(:)); + setCount(k, 1) = 0; + end + setCount(k) = setCount(k) + 1; + end +end + +% Merge the enzyme sets that overlap enough +nSets = numel(setIdx); +parent = 1:nSets; +byEnzyme = cell(numel(enzymes), 1); +for i = 1:nSets + for e = setIdx{i}' + byEnzyme{e}(end+1) = i; + end +end +for e = 1:numel(byEnzyme) + indices = byEnzyme{e}; + for x = 1:numel(indices) + for y = x+1:numel(indices) + i = indices(x); + j = indices(y); + common = numel(intersect(setIdx{i}, setIdx{j})); + if common >= minShared && common >= minOverlap * max(numel(setIdx{i}), numel(setIdx{j})) + parent(findRoot(parent, i)) = findRoot(parent, j); + end + end + end +end +roots = zeros(1, nSets); +for i = 1:nSets + roots(i) = findRoot(parent, i); +end + +% One complex per group of sets, numbered by the first set of the group +[groupRoots, ~, groupOf] = unique(roots, 'stable'); +nComplex = numel(groupRoots); +complexes = cell(nComplex, 1); +weight = cell(nComplex, 1); +for g = 1:nComplex + idx = zeros(0, 1); + copies = zeros(0, 1); + w = zeros(numel(enzymes), 1); + for i = find(groupOf(:)' == g) + w(setIdx{i}) = w(setIdx{i}) + setCount(i); + for n = 1:numel(setIdx{i}) + e = setIdx{i}(n); + pos = find(idx == e, 1); + if isempty(pos) + idx(end+1, 1) = e; %#ok + copies(end+1, 1) = setCopies{i}(n); %#ok + elseif copies(pos) ~= setCopies{i}(n) + warning('getComplexFamilies:copiesDisagree', ... + '%s has %g copies in one variant and %g in another; the first is used.', ... + enzymes{e}, copies(pos), setCopies{i}(n)) + end + end + end + complexes{g} = struct('enzymes', {enzymes(idx)}, 'copies', copies, 'idx', idx); + weight{g} = w; +end + +% A subunit belongs to the complex with most reactions when it is in two +owner = zeros(numel(enzymes), 1); +inComplex = false(numel(enzymes), 1); +for e = 1:numel(enzymes) + hits = []; + for g = 1:nComplex + if weight{g}(e) > 0 + hits(end+1, :) = [weight{g}(e), g]; %#ok + end + end + if isempty(hits) + continue + end + inComplex(e) = true; + if alone(e) || size(hits, 1) > 2 + continue + end + if size(hits, 1) == 1 + owner(e) = hits(1, 2); + elseif hits(1, 1) ~= hits(2, 1) + [~, best] = max(hits(:, 1)); + owner(e) = hits(best, 2); + end +end +for g = 1:nComplex + complexes{g}.owned = owner(complexes{g}.idx) == g; + complexes{g} = rmfield(complexes{g}, 'idx'); +end + +families.complexes = complexes; +families.owner = owner; +families.shared = enzymes(inComplex & owner == 0); +end + +function r = findRoot(parent, i) +% Root of the union-find tree that contains i +r = i; +while parent(r) ~= r + r = parent(r); +end +end diff --git a/src/limit_proteins/levelComplexConcs.m b/src/limit_proteins/levelComplexConcs.m new file mode 100644 index 00000000..9b927996 --- /dev/null +++ b/src/limit_proteins/levelComplexConcs.m @@ -0,0 +1,132 @@ +function [model, result] = levelComplexConcs(model, varargin) +% levelComplexConcs Give every enzyme complex one abundance. +% +% Absolute proteomics recovers the subunits of a complex unevenly, so after +% fillEnzConcs the subunits of one complex can disagree several fold about +% how much of it there is. Capped one by one, the complex can then carry no +% more flux than its least abundant subunit allows. This function groups +% the enzymes of every multi-subunit reaction into complexes (see +% getComplexFamilies) and gives each complex one abundance, scaled to every +% subunit by its stoichiometry. +% +% 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 +% that are shared between complexes, or that also catalyse a reaction on +% their own, are left as measured. +% +% Parameters +% ---------- +% model : struct +% an ecModel in GECKO 3 format (with ecModel.ec structure), with +% ec.concs filled. +% +% Name-Value Arguments +% -------------------- +% statistic : char +% the summary over the measured subunits, '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). +% families : struct +% the complexes of the model, as returned by getComplexFamilies +% (default: calculated from the model). +% +% Returns +% ------- +% model : struct +% ecModel where model.ec.concs has one abundance per complex. +% result : struct +% what was changed, with the fields: +% +% - adjusted : cell array with the measured subunits whose +% concentration was replaced by that of the complex. +% - added : cell array with the unmeasured subunits that were given the +% concentration of the complex. +% - shared : cell array with the measured subunits that were left as +% 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. +% +% Examples +% -------- +% % optional arguments may be given positionally or as name-value pairs: +% model = levelComplexConcs(model); +% model = levelComplexConcs(model, 'statistic', 'mean'); +% [model, result] = levelComplexConcs(model, 'fillMissing', false); +% +% See also +% -------- +% getComplexFamilies, fillEnzConcs, constrainEnzConcs, flexibilizeEnzConcs + +p = parseGECKOargs(varargin, { ... + 'statistic', []; ... + 'fillMissing', []; ... + 'families', []}); +statistic = p.statistic; +fillMissing = p.fillMissing; +families = p.families; +if isempty(statistic); statistic = 'median'; end +if isempty(fillMissing); fillMissing = true; end +statistic = char(statistic); +if ~ismember(statistic, {'median','mean'}) + error('levelComplexConcs:statistic', ... + 'statistic must be ''mean'' or ''median'', not ''%s''.', statistic) +end +summarise = str2func(statistic); +if isempty(families) + families = getComplexFamilies(model); +end + +if ~isfield(model.ec, 'concs') || numel(model.ec.concs) ~= numel(model.ec.enzymes) + error('levelComplexConcs:noConcs', ... + 'model.ec.concs has no entry per enzyme. Run fillEnzConcs first.') +end +enzymes = model.ec.enzymes; +mw =model.ec.mw(:); +concs = model.ec.concs(:); +molar = concs ./ mw; + +% Measured amounts by enzyme: finite entries, in order of the enzymes +values = molar; +isKnown = isfinite(molar); +isLevel = isKnown & molar > 0; +result.adjusted = cell(0, 1); +result.added = cell(0, 1); +result.shared = sort(enzymes(ismember(enzymes, families.shared) & isKnown)); +result.skipped = cell(0, 1); + +for c = 1:numel(families.complexes) + cplx = families.complexes{c}; + [~, pos] = ismember(cplx.enzymes, enzymes); + owned = find(cplx.owned); + ownedPos = pos(owned); + measured = isLevel(ownedPos); + if ~any(measured) + result.skipped{end+1, 1} = sort(cplx.enzymes(:))'; %#ok + continue + end + perCopy = summarise(values(ownedPos(measured)) ./ cplx.copies(owned(measured))); + for n = 1:numel(owned) + e = ownedPos(n); + level = perCopy * cplx.copies(owned(n)); + if isLevel(e) + if abs(values(e) - level) > 1e-8 + 1e-5 * abs(level) + result.adjusted{end+1, 1} = enzymes{e}; %#ok + end + values(e) = level; + elseif fillMissing + values(e) = level; + isKnown(e) = true; + result.added{end+1, 1} = enzymes{e}; %#ok + end + end +end + +model.ec.concs(isKnown) = values(isKnown) .* mw(isKnown); +end diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index bf792183..253402a0 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -65,6 +65,9 @@ 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/getComplexFamilies.m tc0262, tc0263, tc0264, tc0265, tc0266 +limit_proteins/levelComplexConcs.m tc0267, tc0268, tc0269, tc0270, tc0271, tc0272 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/levelComplexTests.m b/test/unit_tests/levelComplexTests.m new file mode 100644 index 00000000..63905d41 --- /dev/null +++ b/test/unit_tests/levelComplexTests.m @@ -0,0 +1,297 @@ +%run this test case with the command +%results = runtests('levelComplexTests.m') + +function tests = levelComplexTests + tests = functiontests(localfunctions); +end + +function testComplexFamiliesGroupIsozymesAndFlagSharedSubunits_tc0262(testCase) + % B is in an isozyme variant with C, so A, B and C are one complex; S is + % in two complexes of equal weight and belongs to none. + model = ecOnly(nan(9, 1), toyRows()); + families = getComplexFamilies(model); + verifyEqual(testCase, sortedComplexes(families), {{'A','B','C'}, {'D','F','S'}, {'E','G','S'}}) + verifyEqual(testCase, families.shared, {'S'}) + [~, b] = ismember('B', model.ec.enzymes); + [~, s] = ismember('S', model.ec.enzymes); + [~, h] = ismember('H', model.ec.enzymes); + verifyEqual(testCase, families.complexes{families.owner(b)}.enzymes(families.complexes{families.owner(b)}.owned), {'A';'B';'C'}) + verifyEqual(testCase, families.owner(s), 0) + verifyEqual(testCase, families.owner(h), 0) +end + +function testComplexFamiliesEnzymeThatWorksAloneIsShared_tc0263(testCase) + rows = toyRows(); + rows(end+1, :) = {'R5', struct('A', 1)}; + families = getComplexFamilies(ecOnly(nan(9, 1), rows)); + verifyTrue(testCase, ismember('A', families.shared)) +end + +function testComplexFamiliesCarrierDoesNotBridgeComplexes_tc0264(testCase) + % B, C, D take carrier S or carrier H (two variants); S also pairs with + % E, like cytochrome c with a dehydrogenase. A single shared protein does + % not merge the two. + rows = {'R1_EXP_1', struct('B',1,'C',1,'D',1,'S',1); ... + 'R1_EXP_2', struct('B',1,'C',1,'D',1,'H',1); ... + 'R2', struct('E',1,'S',1)}; + families = getComplexFamilies(ecOnly(nan(9, 1), rows)); + verifyEqual(testCase, sortedComplexes(families), {{'B','C','D','H','S'}, {'E','S'}}) + verifyTrue(testCase, ismember('S', families.shared)) +end + +function testComplexFamiliesSubunitInTwoComplexesGoesToTheHeavier_tc0265(testCase) + % A, B, C, D is a complex on three reactions; A, B with E, F, G, H on + % one (an assembly step): A and B belong to the first complex. + big = struct('A',1,'B',1,'C',1,'D',1); + rows = {'R1', big; 'R2', big; 'R3', big; ... + 'R4', struct('A',1,'B',1,'E',1,'F',1,'G',1,'H',1)}; + model = ecOnly(nan(9, 1), rows); + families = getComplexFamilies(model); + verifyEmpty(testCase, families.shared) + verifyEqual(testCase, sort(ownedBy(families, model, 'A')), {'A';'B';'C';'D'}) + verifyEqual(testCase, sort(ownedBy(families, model, 'E')), {'E';'F';'G';'H'}) + % with equal weights the two subunits are shared + rows = {'R1', big; 'R2', struct('A',1,'B',1,'E',1,'F',1,'G',1,'H',1)}; + families = getComplexFamilies(ecOnly(nan(9, 1), rows)); + verifyTrue(testCase, all(ismember({'A','B'}, families.shared))) +end + +function testComplexFamiliesWarnWhenVariantsDisagreeOnCopies_tc0266(testCase) + rows = {'R1', struct('A',1,'B',2); 'R2', struct('A',1,'B',3,'C',1)}; + model = ecOnly(nan(9, 1), rows); + verifyWarning(testCase, @() getComplexFamilies(model), 'getComplexFamilies:copiesDisagree') + warning('off', 'getComplexFamilies:copiesDisagree'); + families = getComplexFamilies(model); + warning('on', 'getComplexFamilies:copiesDisagree'); + [~, b] = ismember('B', families.complexes{1}.enzymes); + verifyEqual(testCase, families.complexes{1}.copies(b), 2) +end + +function testLevelComplexConcsMedianScalesByCopies_tc0267(testCase) + % Molar amount per copy: A 1, B 8 / 2 = 4, C 2; the median 2 gives + % A = 2 * 1 copy, B = 2 * 2 copies, C = 2 * 1 copy. + concs = nan(9, 1); + concs(1) = 100; concs(2) = 1600; concs(3) = 100; + model = ecOnly(concs, toyRows()); + [model, result] = levelComplexConcs(model); + verifyEqual(testCase, model.ec.concs(1:3), [200; 800; 100], 'RelTol', 1e-12) + verifyTrue(testCase, isnan(model.ec.concs(9))) + verifyEqual(testCase, sort(result.adjusted), {'A';'B'}) + verifyEmpty(testCase, result.added) +end + +function testLevelComplexConcsMeanIsAvailable_tc0268(testCase) + % Molar amount per copy: A 1, B 4 / 2 = 2, C 6; the mean is 3. + concs = nan(9, 1); + concs(1) = 100; concs(2) = 800; concs(3) = 300; + model = levelComplexConcs(ecOnly(concs, toyRows()), 'statistic', 'mean'); + verifyEqual(testCase, model.ec.concs(1:3), [300; 1200; 150], 'RelTol', 1e-12) + model = levelComplexConcs(ecOnly(concs, toyRows()), 'median'); + verifyEqual(testCase, model.ec.concs(1:3), [200; 800; 100], 'RelTol', 1e-12) +end + +function testLevelComplexConcsSharedSubunitStaysAndMissingOneIsFilled_tc0269(testCase) + concs = nan(9, 1); + concs(4) = 100; % D + concs(6) = 900; % S, shared between two complexes + [model, result] = levelComplexConcs(ecOnly(concs, toyRows())); + verifyEqual(testCase, model.ec.concs([4 5 6]), [100; 100; 900]) + verifyEqual(testCase, result.shared, {'S'}) + verifyEqual(testCase, result.added, {'F'}) + % filling can be switched off + model = levelComplexConcs(ecOnly(concs, toyRows()), 'fillMissing', false); + verifyTrue(testCase, isnan(model.ec.concs(5))) + verifyEqual(testCase, model.ec.concs(6), 900) +end + +function testLevelComplexConcsFillsSubunitOfIsozymeVariant_tc0270(testCase) + % C is in the second variant and is not measured: it gets the + % median per copy (A 1, B 4 / 2 = 2: 1.5) at its molecular weight. + concs = nan(9, 1); + concs(1) = 100; concs(2) = 800; + [model, result] = levelComplexConcs(ecOnly(concs, toyRows())); + verifyEqual(testCase, model.ec.concs(1:3), [150; 600; 75], 'RelTol', 1e-12) + verifyEqual(testCase, result.added, {'C'}) +end + +function testLevelComplexConcsReportsComplexWithoutMeasurement_tc0271(testCase) + concs = nan(9, 1); + concs(1) = 100; % complex A, B, C is measured + [model, result] = levelComplexConcs(ecOnly(concs, toyRows())); + verifyEqual(testCase, result.skipped, {{'D','F','S'}; {'E','G','S'}}) + verifyTrue(testCase, all(isnan(model.ec.concs(4:9)))) +end + +function testLevelComplexConcsRejectsUnknownStatistic_tc0272(testCase) + model = ecOnly(nan(9, 1), toyRows()); + verifyError(testCase, @() levelComplexConcs(model, 'statistic', 'max'), 'levelComplexConcs:statistic') +end + +function testFlexibilizePerComplexRaisesSubunitsByTheSameFactor_tc0273(testCase) + adapter = getAdapter(); + % Without perComplex the limiting subunit alone is raised + [model, flexEnz] = flexibilizeEnzConcs(complexToy([1 3], [1 3]), 2, [], [], adapter, false); + verifyEqual(testCase, flexEnz.uniprotIDs, {'X1'}) + verifyEqual(testCase, usageUb(model, 'X2'), 3) + % With perComplex both are raised by one factor + [model, flexEnz] = flexibilizeEnzConcs(complexToy([1 3], [1 3]), 'expGrowth', 2, ... + 'modelAdapter', adapter, 'verbose', false, 'perComplex', true); + verifyEqual(testCase, sort(flexEnz.uniprotIDs), {'X1';'X2'}) + [~, order] = ismember({'X1';'X2'}, flexEnz.uniprotIDs); + verifyEqual(testCase, flexEnz.flexConcs(order), [1.995; 5.985], 'RelTol', 1e-6) + verifyEqual(testCase, usageUb(model, 'X1') / 1, usageUb(model, 'X2') / 3, 'RelTol', 1e-9) + sol = solveLP(model); + verifyEqual(testCase, sol.f, 2, 'RelTol', 5e-3) + % a cap is never lowered below its measurement + verifyGreaterThanOrEqual(testCase, usageUb(model, 'X1'), 1) + verifyGreaterThanOrEqual(testCase, usageUb(model, 'X2'), 3) + verifyEqual(testCase, flexEnz.ratioIncr, flexEnz.flexConcs ./ flexEnz.oldConcs) +end + +function testFlexibilizePerComplexKeepsLevelledComplexLevelled_tc0274(testCase) + % Both subunits capped at 1 bind together: no single subunit shows a gain + adapter = getAdapter(); + [model, flexEnz] = flexibilizeEnzConcs(complexToy([1 1], [1 1]), 'expGrowth', 2, ... + 'modelAdapter', adapter, 'verbose', false, 'perComplex', true); + verifyEqual(testCase, usageUb(model, 'X1'), usageUb(model, 'X2'), 'RelTol', 1e-9) + verifyEqual(testCase, sort(flexEnz.uniprotIDs), {'X1';'X2'}) + sol = solveLP(model); + verifyEqual(testCase, sol.f, 2, 'RelTol', 5e-3) +end + +function testFlexibilizePerComplexFalseIsTheDefault_tc0275(testCase) + adapter = getAdapter(); + [m0, f0] = flexibilizeEnzConcs(complexToy([1 3], [1 3]), 'expGrowth', 2, 'modelAdapter', adapter, 'verbose', false); + [m1, f1] = flexibilizeEnzConcs(complexToy([1 3], [1 3]), 'expGrowth', 2, 'modelAdapter', adapter, 'verbose', false, 'perComplex', false); + verifyEqual(testCase, m1, m0) + verifyEqual(testCase, f1, f0) + verifyEqual(testCase, f0.uniprotIDs, {'X1'}) +end + +function testFlexibilizePerComplexDoesNotPickComplexForPoolPrice_tc0276(testCase) + % A -> B by the complex X1, X2, X3 and B -> P by the single enzyme E, + % with a protein pool that binds. At flux 1 every subunit uses 80% of a + % 1.25 cap, E uses all of its cap of 1. Of the enzymes, E alone limits. + adapter = getAdapter(); + [model, flexEnz] = flexibilizeEnzConcs(seriesToy(), 'expGrowth', 1.1, 'modelAdapter', adapter, ... + 'verbose', false, 'perComplex', true); + verifyEqual(testCase, flexEnz.uniprotIDs(:), {'prot_pool';'E'}) + verifyEqual(testCase, usageUb(model, 'E'), 1.09725, 'RelTol', 1e-4) + for e = {'X1','X2','X3'} + verifyEqual(testCase, usageUb(model, e{1}), 1.25, 'RelTol', 1e-9) + end + sol = solveLP(model); + verifyGreaterThanOrEqual(testCase, sol.f, 1.1 * 0.9975 - 1e-9) +end + +% ---- helpers + +function adapter = getAdapter() + adapter = ModelAdapterManager.getAdapter(fullfile(findGECKOroot,'test','unit_tests','ecTestGEM','TestGEMAdapter.m')); +end + +function rows = toyRows() + rows = {'R1_EXP_1', struct('A',1,'B',2); ... + 'R1_EXP_2', struct('A',1,'B',2,'C',1); ... % isozyme variant: same complex + 'R2', struct('D',1,'F',1,'S',1); ... + 'R3', struct('E',1,'G',1,'S',1); ... % S is shared with the complex of R2 + 'R4', struct('H',1)}; % a single enzyme: not a complex +end + +function model = ecOnly(concs, rows) + % A model that carries nothing but the ec data: rows is {rxn, struct(enzyme, copies)} + enzymes = {'A';'B';'C';'D';'F';'S';'E';'G';'H'}; + mw = [100; 200; 50; 100; 100; 100; 100; 100; 100]; + mat = zeros(size(rows, 1), numel(enzymes)); + for r = 1:size(rows, 1) + names = fieldnames(rows{r, 2}); + for n = 1:numel(names) + mat(r, strcmp(enzymes, names{n})) = rows{r, 2}.(names{n}); + end + end + model.ec.rxns = rows(:, 1); + model.ec.genes = strcat('g_', enzymes); + model.ec.enzymes = enzymes; + model.ec.mw = mw; + model.ec.concs = concs(:); + model.ec.rxnEnzMat = sparse(mat); +end + +function sets = sortedComplexes(families) + sets = cellfun(@(c) sort(c.enzymes(:))', families.complexes, 'UniformOutput', false); + [~, order] = sort(cellfun(@(c) strjoin(c, ''), sets, 'UniformOutput', false)); + sets = sets(order)'; +end + +function names = ownedBy(families, model, enzyme) + [~, e] = ismember(enzyme, model.ec.enzymes); + c = families.complexes{families.owner(e)}; + names = c.enzymes(c.owned); +end + +function ub = usageUb(model, enzyme) + ub = model.ub(strcmp(model.rxns, ['usage_prot_' enzyme])); +end + +function model = complexToy(concs, ubs) + % A -> B is catalysed by the complex X1 + X2; growth (R4) equals the flux. + model = buildLP({'A_e';'A_c';'B_c';'prot_pool';'prot_X1';'prot_X2'}, ... + {'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)}); + 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 + +function model = seriesToy() + % A -> B by the complex X1, X2, X3, B -> P by the enzyme E; the pool + % holds exactly 4 mg. + model = buildLP({'A';'B';'P';'prot_pool';'prot_X1';'prot_X2';'prot_X3';'prot_E'}, ... + {'EX_A', {'A',1}, 0, 1000; ... + 'R1', {'A',-1,'B',1,'prot_X1',-1,'prot_X2',-1,'prot_X3',-1}, 0, 1000; ... + 'R2', {'B',-1,'P',1,'prot_E',-1}, 0, 1000; ... + 'R4', {'P',-1}, 0, 1000; ... + 'prot_pool_exchange', {'prot_pool',1}, 0, 4; ... + 'usage_prot_X1', {'prot_pool',-1,'prot_X1',1}, 0, 1.25; ... + 'usage_prot_X2', {'prot_pool',-1,'prot_X2',1}, 0, 1.25; ... + 'usage_prot_X3', {'prot_pool',-1,'prot_X3',1}, 0, 1.25; ... + 'usage_prot_E', {'prot_pool',-1,'prot_E',1}, 0, 1}); + model.ec.enzymes = {'X1';'X2';'X3';'E'}; + model.ec.rxns = {'R1';'R2'}; + model.ec.genes = {'g1';'g2';'g3';'gE'}; + model.ec.mw = [100; 100; 100; 100]; + model.ec.concs = [1.25; 1.25; 1.25; 1]; + model.ec.rxnEnzMat = sparse([1 1 1 0; 0 0 0 1]); +end + +function model = buildLP(mets, rxnSpec) + % rxnSpec rows: id, {metabolite, coefficient, ...}, lb, ub. The reaction + % named R4 is the objective, which is the biomass reaction of the adapter. + nR = size(rxnSpec, 1); + S = zeros(numel(mets), nR); + for r = 1:nR + entries = rxnSpec{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 = rxnSpec(:, 1); + model.rxnNames = rxnSpec(:, 1); + model.S = sparse(S); + model.lb = cell2mat(rxnSpec(:, 3)); + model.ub = cell2mat(rxnSpec(:, 4)); + model.c = double(strcmp(model.rxns, 'R4')); + model.b = zeros(numel(mets), 1); + model.rev = double(model.lb < 0); +end