Edge cases of GECKO agree with geckopy

src/change_model/addNewRxnsToEC.m

@@ -50,7 +50,9 @@
 %    Note that the reversibility defined in the equation will be used to
 %    split to construct irreversible reactions (add _REV) and the rules
 %    defined will be used to expand the model (add _EXP_n).
-% 2. After adding the new reactions, setKcatForReactions or
+% 2. When all the genes of the new reactions already have enzymes, no
+%    newEnzymes are needed.
+% 3. After adding the new reactions, setKcatForReactions or
 %    applyCustomKcats should be run to add the kcat to the reactions, and
 %    subsequently applyKcatConstraints.
 %
@@ -78,8 +80,9 @@
 end
 params = modelAdapter.getParameters();
 
+enzAdded = {};
 if isempty(newEnzymes)
-    newEnzymes = [];
+    newEnzymes = struct('enzymes', {{}}, 'genes', {{}}, 'mw', []);
 end
 
 % Validate if GECKO version 3
@@ -91,7 +94,7 @@
     error('addNewRxnsToEC only works on non-light ecModels')
 else
     % Check if the new proteins are not already present in the model
-    if ~isempty(newEnzymes)
+    if ~isempty(newEnzymes.enzymes)
         [toRemove,~] = ismember(newEnzymes.enzymes, model.ec.enzymes);
 
         % Remove from the list those enzymes that already present
@@ -213,19 +216,21 @@
     model = addRxns(model, newRxns, 3, '', true, true);
 
     % 4. Add protein usage reactions.
-    usageRxns.rxns            = strcat('usage_',proteinMets.mets);
-    usageRxns.rxnNames        = usageRxns.rxns;
-    usageRxns.mets            = cell(numel(usageRxns.rxns),1);
-    usageRxns.stoichCoeffs    = cell(numel(usageRxns.rxns),1);
-    for i=1:numel(usageRxns.mets)
-        usageRxns.mets{i}         = {proteinMets.mets{i}, 'prot_pool'};
-        usageRxns.stoichCoeffs{i} = [1,-1];
+    if ~isempty(newEnzymes.enzymes)
+        usageRxns.rxns            = strcat('usage_',proteinMets.mets);
+        usageRxns.rxnNames        = usageRxns.rxns;
+        usageRxns.mets            = cell(numel(usageRxns.rxns),1);
+        usageRxns.stoichCoeffs    = cell(numel(usageRxns.rxns),1);
+        for i=1:numel(usageRxns.mets)
+            usageRxns.mets{i}         = {proteinMets.mets{i}, 'prot_pool'};
+            usageRxns.stoichCoeffs{i} = [1,-1];
+        end
+        usageRxns.lb              = zeros(numel(usageRxns.rxns),1);
+        usageRxns.ub              = zeros(numel(usageRxns.rxns),1) + 1000;
+        usageRxns.rev             = ones(numel(usageRxns.rxns),1);
+        usageRxns.grRules         = newEnzymes.genes;
+        model = addRxns(model,usageRxns);
     end
-    usageRxns.lb              = zeros(numel(usageRxns.rxns),1);
-    usageRxns.ub              = zeros(numel(usageRxns.rxns),1) + 1000;
-    usageRxns.rev             = ones(numel(usageRxns.rxns),1);
-    usageRxns.grRules         = newEnzymes.genes;
-    model = addRxns(model,usageRxns);
 
     % 5. Update .ec field with enzyme data
     model.ec.enzymes    = [model.ec.enzymes; newEnzymes.enzymes(:)];

src/change_model/applyComplexData.m

@@ -38,7 +38,9 @@
 % Notes
 % -----
 % The subunit stoichiometry is assigned per protein ID: the proteins of a
-% complex need not be listed in the order of model.ec.enzymes.
+% complex need not be listed in the order of model.ec.enzymes. A subunit with
+% stoichiometry 0 has no stated stoichiometry; it counts as a subunit with
+% stoichiometry 1.
 %
 % foundComplex and proposedComplex share the following columns: rxn
 % (reaction identifier), complexID and name (Complex Portal identifiers),
@@ -113,9 +115,8 @@
     protMatIdx(~isOld)  = lastProt+1 : lastProt+numel(newProts);
     lastProt = lastProt + numel(newProts);
     % Add subunit numbers
-    if all(complexData(i).stochiometry == 0)
-        complexData(i).stochiometry(:) = 1;
-    end
+    % A subunit without a stated stoichiometry (0) counts as one subunit
+    complexData(i).stochiometry(complexData(i).stochiometry == 0) = 1;
     complexMatrix(i,protMatIdx) = complexData(i).stochiometry;
 end
 

src/change_model/applyKcatConstraints.m

@@ -3,9 +3,8 @@
 %
 % Existing enzyme constraints are first removed (unless updateRxns is
 % provided), and new constraints are defined based on the content of
-% model.ec.kcat. When every kcat is 0, the existing constraints of a full
-% ecModel are removed and no new ones are written; a light ecModel is
-% returned unchanged.
+% model.ec.kcat. When every kcat is 0, the existing constraints are removed
+% and no new ones are written.
 %
 % Parameters
 % ----------
@@ -18,8 +17,13 @@
 %     if not all enzyme constraints should be updated, this can be given as
 %     either a logical vector of length model.ec.rxns, a vector of
 %     model.ec.rxns indices, or a (cell array of) string(s) with
-%     model.ec.rxns identifiers. For light models, these reactions should
-%     match model.rxns.
+%     model.ec.rxns identifiers. For light models, a (cell array of) string(s)
+%     may hold reaction identifiers of model.rxns or of model.ec.rxns (with
+%     the ###_ prefix); a reaction is then updated as a whole, with all its
+%     isozymes. An identifier that is in both or in neither raises the
+%     error applyKcatConstraints:ambiguousReaction or
+%     applyKcatConstraints:unknownReaction. A logical vector or indices
+%     refer to model.rxns for light models.
 %
 % Returns
 % -------
@@ -56,6 +60,23 @@
     updateRxns = updateRxnsLog;
 elseif iscellstr(updateRxns) || ischar(updateRxns) || isstring(updateRxns)
     updateRxnsIds = convertCharArray(updateRxns);
+    updateRxnsIds = updateRxnsIds(:);
+    if model.ec.geckoLight
+        isCobra = ismember(updateRxnsIds,model.rxns);
+        isEc    = ismember(updateRxnsIds,model.ec.rxns);
+        if any(isCobra & isEc)
+            error('applyKcatConstraints:ambiguousReaction', ...
+                'The identifiers %s are in both model.rxns and model.ec.rxns.', ...
+                strjoin(updateRxnsIds(isCobra & isEc),', '))
+        end
+        if any(~isCobra & ~isEc)
+            error('applyKcatConstraints:unknownReaction', ...
+                'The identifiers %s are in neither model.rxns nor model.ec.rxns.', ...
+                strjoin(updateRxnsIds(~isCobra & ~isEc),', '))
+        end
+        updateRxnsIds = [updateRxnsIds(isCobra); extractAfter(updateRxnsIds(isEc),4)];
+        updateRxnsIds = updateRxnsIds(:);
+    end
     updateRxns = ismember(rxns,updateRxnsIds);
 end
 
@@ -68,12 +89,8 @@
            ' not a valid GECKO3 ecModel. First run makeEcModel(model).'])
 end
 if all(model.ec.kcat==0)
-    if model.ec.geckoLight
-        printOrange('WARNING: No kcat values are provided in model.ec.kcat, model remains unchanged.\n');
-        return
-    end
-    % In a full model the existing enzyme constraints of the selected
-    % reactions are still cleared below, and no new ones are written.
+    % The existing enzyme constraints of the selected reactions are still
+    % cleared below, and no new ones are written.
     printOrange('WARNING: No kcat values are provided in model.ec.kcat, existing enzyme constraints are cleared and no new ones are written.\n');
 end
 

src/change_model/findMetSmiles.m

@@ -3,7 +3,10 @@
 %
 % Matches are cached in tutorials/***/data/smilesDB.tsv, which is checked
 % first on subsequent runs so repeated queries are avoided. If the model
-% already has a metSmiles field, non-empty entries are left unchanged.
+% already has a metSmiles field, non-empty entries are left unchanged. When
+% the cache holds a name more than once, the last row is used. A metabolite
+% without a name is skipped with a warning (findMetSmiles:emptyName) and
+% gets no SMILES.
 %
 % Parameters
 % ----------
@@ -39,6 +42,10 @@
 end
 params = modelAdapter.params;
 
+nEmpty = sum(cellfun(@isempty,model.metNames));
+if nEmpty > 0
+    warning('findMetSmiles:emptyName','%d metabolite(s) have no name and get no SMILES.',nEmpty)
+end
 [uniqueNames, ~, uniqueIdx] = unique(regexprep(model.metNames,'^prot_.*',''));
 uniqueSmiles(1:numel(uniqueNames),1) = {''};
 metMatch = false(length(uniqueNames),1);
@@ -49,8 +56,9 @@
     fID = fopen(smilesDBfile,'r');
     raw = textscan(fID,'%s %s','Delimiter','\t','HeaderLines',0);
     fclose(fID);
-    smilesDB.names = raw{1};
-    smilesDB.smile = raw{2};
+    %The last row of a name that occurs more than once is used
+    smilesDB.names = flipud(raw{1});
+    smilesDB.smile = flipud(raw{2});
     [metMatch, metIdx] = ismember(uniqueNames,smilesDB.names);
     uniqueSmiles(metMatch) = smilesDB.smile(metIdx(metMatch));
     if verbose; fprintf('done.\n'); end

src/enzyme_usage/getEnzymeUsage.m

@@ -39,7 +39,9 @@
 % - capUsage : vector of enzyme capacity usages.
 % - absUsage : vector of absolute enzyme usages.
 % - UB : vector of enzyme exchange reaction upper bounds.
-% - protID : string array of matching protein IDs.
+% - protID : string array of matching protein IDs. An enzyme without a
+%   usage_prot_<id> reaction in the model is left out, with the warning
+%   getEnzymeUsage:missingUsageReaction.
 % - fluxes : vector of fluxes, copy of input fluxes.
 %
 % Examples
@@ -60,6 +62,14 @@
 end
 usageData.protID      = ecModel.ec.enzymes;
 [~,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);
+    rxnIdx = rxnIdx(~noUsageRxn);
+end
 
 usageData.UB          = ecModel.ub(rxnIdx);
 usageData.absUsage    = abs(fluxes(rxnIdx));

src/enzyme_usage/reportEnzymeUsage.m

@@ -20,7 +20,7 @@
 %     how many top enzymes by absolute usage to list (default 10). Use 0
 %     or Inf to list all enzymes instead. Enzymes whose reactions carry no
 %     flux are never listed here, as they do not contribute to the
-%     solution.
+%     solution. A negative value raises reportEnzymeUsage:invalidTopAbsUsage.
 %
 % Returns
 % -------
@@ -54,6 +54,9 @@
 % values it used to reach the indexing below unclamped: usageData.protID(topUse(1:N))
 % with N > numel(topUse) is an out-of-bounds index, so any model with fewer than 10
 % enzymes (ecTestGEM has 5) crashed on the default call.
+if ~(isnumeric(topAbsUsage) && isscalar(topAbsUsage) && topAbsUsage >= 0)
+    error('reportEnzymeUsage:invalidTopAbsUsage', 'topAbsUsage must be a non-negative number.')
+end
 if topAbsUsage == 0 || isinf(topAbsUsage) || topAbsUsage > numel(usageData.protID)
     topAbsUsage = numel(usageData.protID);
 end
@@ -84,8 +87,8 @@
     if isscalar(find(carriedFlux))
         highUsage.protID(end+1,1)      = highEnzyme(i);
         highUsage.geneID(end+1,1)      = geneIDs(i);
-        highUsage.absUsage(end+1,1)    = usageData.absUsage(enzIdx(i));
-        highUsage.capUsage(end+1,1)    = usageData.capUsage(enzIdx(i));
+        highUsage.absUsage(end+1,1)    = usageData.absUsage(highUsageProt(i));
+        highUsage.capUsage(end+1,1)    = usageData.capUsage(highUsageProt(i));
         highUsage.kcat(end+1,1)        = kcat(carriedFlux);
         highUsage.source(end+1,1)      = ecModel.ec.source(idx(carriedFlux));
         highUsage.rxnID(end+1,1)       = rxns(carriedFlux);
@@ -95,8 +98,8 @@
         % Add one entry for combined usage
         highUsage.protID(end+1,1)      = highEnzyme(i);
         highUsage.geneID(end+1,1)      = geneIDs(i);
-        highUsage.absUsage(end+1,1)    = usageData.absUsage(enzIdx(i));
-        highUsage.capUsage(end+1,1)    = usageData.capUsage(enzIdx(i));
+        highUsage.absUsage(end+1,1)    = usageData.absUsage(highUsageProt(i));
+        highUsage.capUsage(end+1,1)    = usageData.capUsage(highUsageProt(i));
         highUsage.kcat(end+1,1)        = nan;
         highUsage.source{end+1,1}      = '===';
         highUsage.rxnID{end+1,1}       = '===';
@@ -109,7 +112,7 @@
         enzMet = strcat('prot_',highEnzyme{i});
         [~, enzEcIdx] = ismember(enzMet,ecModel.mets);
         indAbsUse = full(transpose(-ecModel.S(enzEcIdx,rIdx)).*enzFlux);
-        indCapUse = (indAbsUse /sum(indAbsUse)) * usageData.capUsage(enzIdx(i));
+        indCapUse = (indAbsUse /sum(indAbsUse)) * usageData.capUsage(highUsageProt(i));
 
         rxnNumber = length(rIdx);
         highUsage.protID(end+1:end+rxnNumber,1)      = highEnzyme(i);

src/flux_analysis/ecFSEOF.m

@@ -194,8 +194,10 @@
 corr_rxns   = corr_rxns(targets);
 target_type = target_type(targets);
 
-% Order from highest to lowest score, the product of slope and correlation
-[~,order]   = sort(slope_rxns .* corr_rxns,'descend');
+% Order from highest to lowest score, the product of slope and correlation,
+% rounded to 9 significant digits so that scores equal up to solver noise keep
+% the order of the model
+[~,order]   = sort(round(slope_rxns .* corr_rxns,9,'significant'),'descend');
 target_rxns = target_rxns(order);
 v_matrix    = v_matrix(order,:);
 rxnGeneM    = rxnGeneM(order,:);

src/gather_kcats/assignKcatValues.m

@@ -33,6 +33,11 @@
 %     'median'/'mean', where the selected value need not equal any single
 %     input, it is taken from the first entry instead (arbitrary but
 %     deterministic, since an aggregate has no single "winning" source).
+%     The source is written in lower case snake_case, with the wildcard
+%     level and origin of the entry in brackets when it has them
+%     (kcatList.wildcardLvl and kcatList.origin), for example
+%     'brenda (wc=0, origin=1)'. Functions that compare model.ec.source with
+%     a source name use the part before ' (' (see resolveTunableContext).
 % overwrite : char
 %     whether existing kcat values should be overwritten. Options: 'true',
 %     'false', 'ifHigher'. The last option will overwrite only if the new
@@ -77,6 +82,17 @@
 
 % Remove zero kcat values. Only adjusting fields that are used later.
 removeZero                      = kcatList.kcats == 0;
+nEntries                        = numel(removeZero);
+wildcardLvl = nan(nEntries, 1);
+origin      = nan(nEntries, 1);
+if isfield(kcatList,'wildcardLvl') && numel(kcatList.wildcardLvl) == nEntries
+    wildcardLvl = double(kcatList.wildcardLvl(:));
+end
+if isfield(kcatList,'origin') && numel(kcatList.origin) == nEntries
+    origin = double(kcatList.origin(:));
+end
+wildcardLvl(removeZero)         = [];
+origin(removeZero)              = [];
 kcatList.kcats(removeZero)      = [];
 kcatList.rxns(removeZero)       = [];
 if isfield(kcatList,'kcatSource')
@@ -114,7 +130,8 @@
         otherwise
             error('Invalid criteria specified')
     end
-    selectedSource(i)    = kcatList.kcatSource(idxMatch(j));
+    selectedSource{i}    = formatKcatSource(kcatList.kcatSource{idxMatch(j)}, ...
+        wildcardLvl(idxMatch(j)), origin(idxMatch(j)));
 end
 
 % Populate model.ec.kcat
@@ -140,3 +157,17 @@
 end
 rxnIdx = idxInModelUnique;
 end
+
+function text = formatKcatSource(source, wildcardLvl, origin)
+% The provenance text of model.ec.source: the source label in lower case
+% snake_case, followed by the wildcard level and origin of a fuzzy BRENDA
+% match in brackets, e.g. 'brenda (wc=0, origin=1)'.
+text = lower(strtrim(char(source)));
+text = regexprep(text, '[^0-9a-z]+', '_');
+text = regexprep(text, '^_+|_+$', '');
+if ~isnan(wildcardLvl) && ~isnan(origin)
+    text = sprintf('%s (wc=%d, origin=%d)', text, fix(wildcardLvl), fix(origin));
+elseif ~isnan(wildcardLvl)
+    text = sprintf('%s (wc=%d)', text, fix(wildcardLvl));
+end
+end

src/gather_kcats/collectKcats.m

@@ -198,6 +198,7 @@
         continue
     end
     lists = {};
+    usedFiles = {};
     priority = {};
     for t = 1:numel(st.files)
         token = st.files{t};
@@ -209,8 +210,9 @@
         if isempty(list)
             continue
         end
-        if ~any(cellfun(@(l) isequal(l, list), lists))
+        if ~ismember(name, usedFiles)
             lists{end+1} = list; %#ok<AGROW>
+            usedFiles{end+1} = name; %#ok<AGROW>
         end
         if startsWith(token, 'database_')
             priority{end+1} = token; %#ok<AGROW>
@@ -284,10 +286,12 @@
 end
 
 function s = normalizeSource(label)
-% As mergeKcats folds source labels: lower case, other characters as '_'.
+% As mergeKcats folds source labels: lower case, other characters as '_',
+% without a leading okp_ (the tag of OpenKineticsPredictor rows).
 s = lower(strtrim(char(label)));
 s = regexprep(s, '[^0-9a-z]+', '_');
 s = regexprep(s, '^_+|_+$', '');
+s = regexprep(s, '^okp_', '');
 end
 
 function name = sourceFile(token)