Measured zero concentrations and isozyme coverage of the proteomics

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)<maxMisRepl*size(condAbund,2);
+            if zerosAsMeasured
+                present = condAbund>=0;
+            else
+                present = condAbund>0;
+            end
+            remData = sum(present,2)<maxMisRepl*size(condAbund,2);
             condAbund(remData,:) = nan;
         end
         % Filter by RSD
@@ -159,12 +192,19 @@
         condAbund = mean(condAbund,2,'omitnan')+(addStdevs*std(condAbund,0,2,'omitnan'));
         % Filter by minVal
         remData = mean(condAbund,2,'omitnan')<minVal;
+        if zerosAsMeasured
+            remData = remData & condAbund~=0;
+        end
         condAbund(remData) = nan;
         % Remove the lowest cutLowest percent
-        numData  = find(~isnan(condAbund));
-        [~,sortData] = sort(condAbund);
-        lowCutoff = floor(numel(numData)*cutLowest/100);
-        condAbund(sortData(1:lowCutoff)) = nan;
+        if zerosAsMeasured
+            cutData = find(~isnan(condAbund) & condAbund~=0);
+        else
+            cutData = find(~isnan(condAbund));
+        end
+        [~,sortData] = sort(condAbund(cutData));
+        lowCutoff = floor(numel(cutData)*cutLowest/100);
+        condAbund(cutData(sortData(1:lowCutoff))) = nan;
         % Combine conditions
         filtAbund(:,i) = condAbund;
     end

src/limit_proteins/constrainEnzConcs.m

@@ -28,6 +28,12 @@
 % 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.
 %
+% A NaN entry means that the enzyme was not measured, and its usage draws
+% from the protein pool. An entry of 0 means that the enzyme was measured as
+% absent: its usage upper bound is 0, which removes every flux through the
+% enzyme. flexibilizeEnzConcs multiplies the concentration and cannot raise a
+% bound of 0; relaxProteomicsGreedy sets the bound back to its default.
+%
 % Examples
 % --------
 %     % optional arguments may be given positionally or as name-value pairs:

src/limit_proteins/fillEnzConcs.m

@@ -35,6 +35,12 @@
 % A negative abundance of a protein that is an enzyme of the model is
 % contradictory data and raises the error fillEnzConcs:negativeAbundance.
 %
+% A NaN entry in model.ec.concs means that the protein was not measured, and
+% the usage of the enzyme draws from the protein pool. An entry of 0 means
+% that the protein was measured as absent, and constrainEnzConcs limits the
+% usage of the enzyme to 0. loadProtData returns a measured 0 only with the
+% option zerosAsMeasured.
+%
 % To also constrain the model with the content of model.ec.concs, you should
 % run constrainEnzConcs.
 %

src/limit_proteins/flexibilizeEnzConcs.m

@@ -49,7 +49,7 @@
 %     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.
+%     complex is recorded on its first subunit with a concentration above 0.
 %
 % Returns
 % -------
@@ -69,6 +69,15 @@
 %     - frequence : numeric how often the enzyme has been step-wise
 %       flexibilized.
 %
+% Notes
+% -----
+% The usage upper bound of an enzyme is raised by a factor of its
+% concentration, so an enzyme with a concentration of 0 in model.ec.concs
+% stays at 0 and is never flexibilized. If such enzymes exist, the warning
+% flexibilizeEnzConcs:measuredZero lists them: a growth rate that is not
+% reached can then be caused by these enzymes, and flexibilizeEnzConcs
+% reports no limiting enzyme for them.
+%
 % Examples
 % --------
 %     % optional arguments may be given positionally or as name-value pairs:
@@ -143,6 +152,17 @@
 
 % Get those proteins with a concentration defined
 protConcs = find(~isnan(model.ec.concs));
+measuredZero = model.ec.enzymes(model.ec.concs == 0);
+if ~isempty(measuredZero)
+    listed = strjoin(measuredZero(1:min(10, end))', ', ');
+    if numel(measuredZero) > 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);

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_<n> 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_<n> 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

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<AGROW>
             end
             values(e) = level;
-        elseif fillMissing
+        elseif fillMissing && ~isZero(e)
             values(e)  = level;
             isKnown(e) = true;
             result.added{end+1, 1} = enzymes{e}; %#ok<AGROW>

src/limit_proteins/relaxProteomicsGreedy.m

@@ -5,7 +5,9 @@
 % While growth is below minimalGrowth, finds the still-proteomics-constrained
 % enzyme (usage_prot_<id> upper bound below defaultUpperBound) with the
 % largest absolute shadow price on its prot_<id> 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
 % ----------