Fix four bugs found in the help-block review, correct examples and comments

CHANGELOG.md

@@ -21,6 +21,27 @@
 * `gapFillTopological`: without `targets` and without an objective reaction it raises
   `RAVEN:badInput` (was: an error without identifier). Help states the default seeds and
   targets.
+* `checkProduction` with `checkNeededForProduction`: the returned `model` holds excretion
+  reactions only (was: the excretion reaction of the last metabolite in `notProduced` was
+  returned as an uptake reaction). This matters for code that uses the fifth output.
+
+### Curation
+
+* `renameModelGenes` with a file path as `geneTable`: the table is read and the genes are
+  renamed (was: an error `MATLAB:invalidConversion` when the message that names the file
+  was printed).
+
+### Queries
+
+* `getIndexes` with type `'metnames'` and `returnLogical` true: raises `RAVEN:badInput` (was:
+  the error "Unable to use a value of type cell as an index"). A name can match several
+  metabolites, so no single logical vector is defined.
+
+### Analysis
+
+* `findGeneDeletions` with `'dgd'`: `grRatioMuts` has one element per gene pair, as the help
+  states (was: twice as many, the second half zero). This matters for code that indexes
+  `grRatioMuts` by the number of its elements.
 
 ### Context-specific models (ftINIT / tINIT)
 
@@ -57,6 +78,10 @@
   case, `&`, `&&`, `|`, `||`, operators as whole words, so `G1&G2` is one gene id), and that a
   malformed rule raises `RAVEN:badGrRule`.
 
+* `checkRaven`: when the binaries cannot be downloaded, the message names the
+  `RAVEN-<tag>-full.zip` bundle attached to each release for an offline install (was:
+  "the `*-binaries` RAVEN release").
+
 ### Documentation
 
 * Help blocks of 176 files are rewritten in impersonal, present-tense wording with British spelling; the code is unchanged. Statements that differed from the code are corrected, for example `simplifyModel` `deleteDuplicates` (first occurrence kept), `removeBadRxns` `refModel`, `replaceMets` outputs, `getMinNrFluxes` `params`, `checkRaven` `developMode` and the `downloadRavenBinaries` offline bundle. `reverseRxns` has the RAVEN 3 help structure.

checkRaven.m

@@ -257,7 +257,7 @@
         downloadRavenBinaries({'blast+','diamond','hmmer'});
     catch ME
         printOrange(['    Could not download binaries from toolbox-assets: ' ME.message '\n' ...
-            '    (offline? fetch the "*-binaries" RAVEN release for an offline install)\n']);
+            '    (offline? install from the RAVEN-<tag>-full.zip bundle attached to each RAVEN release)\n']);
     end
     % Run each bundled binary with a version/help flag to confirm it
     % executes on this system, rather than running the full test cases.

src/INIT/ftINIT.m

@@ -285,10 +285,10 @@
     lastError = [];
     if ~isempty(proveAbsGap)
         %One solve proven to this fixed absolute gap, replacing the relative-gap
-        %escalation below. The escalation's final near-zero-objective run otherwise
-        %accepts an arbitrary within-gap incumbent -- measured on genome-scale
-        %Human-GEM to cost 2.0-4.0 objective units (raven-toolbox's prove_abs_gap,
-        %raven-gecko-parity#104). Pairs naturally with resolveTies.
+        %escalation below. The final near-zero-objective run of the escalation
+        %accepts an arbitrary incumbent within the gap; on genome-scale Human-GEM
+        %the step 1 objective of the escalation ended 4.0 to 12.0 units below the
+        %proven optimum (raven-toolbox's prove_abs_gap, raven-gecko-parity#104).
         if ~isempty(stp.MILPParams) && isfield(stp.MILPParams{1}, 'TimeLimit')
             params = struct('TimeLimit', stp.MILPParams{1}.TimeLimit);
         else

src/analysis/findGeneDeletions.m

@@ -113,7 +113,7 @@
     [~, I]=ismember(originalGenes(details==1),model.genes);
     genesToModify=nchoosek(I,2);
     genes=geneMapping(genesToModify);
-    grRatioMuts=zeros(1,numel(genes));
+    grRatioMuts=zeros(1,size(genesToModify,1));
     fluxes=sparse(numel(model.rxns),size(genesToModify,1));
     for i=1:size(genesToModify,1)
         tempModel=removeGenes(model,genesToModify(i,:),false,false,false);

src/curation/renameModelGenes.m

@@ -80,14 +80,15 @@
 
     % Load table if a file path was supplied
     if ischar(geneTable) || isstring(geneTable)
-        if ~isfile(geneTable)
+        tableFile = char(geneTable);
+        if ~isfile(tableFile)
             error('renameModelGenes:fileNotFound', ...
-                  'Table file not found: %s', geneTable);
+                  'Table file not found: %s', tableFile);
         end
-        geneTable = readtable(geneTable, 'Delimiter', '\t', 'FileType', 'text', ...
+        geneTable = readtable(tableFile, 'Delimiter', '\t', 'FileType', 'text', ...
                               'TextType', 'char');
         fprintf('Loaded gene table from file: %s  (%d rows)\n', ...
-                geneTable, height(geneTable));
+                tableFile, height(geneTable));
     end
 
     % Validate column names

src/gapfilling/checkProduction.m

@@ -52,11 +52,8 @@
 %     checkNeededForProduction is false.
 % model : struct
 %     model structure with an excretion reaction added for each metabolite
-%     in the compartments of excretionFromCompartments. When
-%     checkNeededForProduction is true, the excretion reaction of the last
-%     metabolite in notProduced is returned as an uptake reaction, because
-%     the loop that flips each reaction into an uptake reaction does not
-%     reset the last one.
+%     in the compartments of excretionFromCompartments. All added reactions
+%     are excretion reactions, also when checkNeededForProduction is true.
 %
 % Examples
 % --------
@@ -104,6 +101,12 @@
         %be produced now
         neededForProductionMat(i,:)=haveFlux(model,10^-5,addedRxns(notProduced));
     end
+    if ~isempty(notProduced)
+        %Reset the last iteration, so that the returned model has excretion
+        %reactions only
+        lastCol=numel(model.rxns)-numel(addedRxns)+notProduced(end);
+        model.S(:,lastCol)=model.S(:,lastCol)*-1;
+    end
     %Calculate the smallest number of metabolites that must be connected to
     %make everything connected and return their names
 

src/gapfilling/fitTasks.m

@@ -65,8 +65,8 @@
 %
 % Examples
 % --------
-%     [outModel, addedRxns]=fitTasks(model,refModel,inputFile,printOutput,...
-%         rxnScores,taskStructure,params);
+%     [outModel, addedRxns, failedTasks]=fitTasks(model,refModel,inputFile,...
+%         'printOutput',false,'rxnScores',rxnScores,'params',params);
 
 p=parseRAVENargs(varargin, {'printOutput',true; ...
     'rxnScores',[]; ...

src/gapfilling/gapReport.m

@@ -69,8 +69,9 @@
 % Examples
 % --------
 %     [noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, ...
-%         minToConnect, neededForProductionMat, connectedFromTemplates, ...
-%         addedFromTemplates] = gapReport(model, templateModels);
+%         minToConnect, neededForProductionMat, canProduceWithoutInput, ...
+%         canConsumeWithoutOutput, connectedFromTemplates, addedFromTemplates] = ...
+%         gapReport(model, 'templateModels', templateModels);
 
 p=parseRAVENargs(varargin, {'templateModels',[]});
 templateModels=p.templateModels;

src/localization/predictLocalization.m

@@ -64,7 +64,8 @@
 % --------
 %     [outModel, geneLocalization, transportStruct, scores, removedRxns] = ...
 %         predictLocalization(model, GSS, defaultCompartment, ...
-%             transportCost, maxTime, plotResults);
+%             'transportCost', transportCost, 'maxTime', maxTime, ...
+%             'plotResults', plotResults);
 
 p=parseRAVENargs(varargin, {'transportCost',[]; 'maxTime',[]; 'plotResults',[]});
 transportCost=p.transportCost;

src/manipulation/closeModel.m

@@ -36,15 +36,9 @@
     col=closedModel.S(:,i);
     %A boundary/exchange/sink/demand reaction has metabolites on only one
     %side, same definition as getExchangeRxns: no substrates or no
-    %products. The previous coefficient-magnitude rule (all nonzero
-    %coefficients summing to 1 in absolute value) missed a scaled
-    %single-metabolite sink such as "2 A =>" (sum 2) outright, and only
-    %matched a multi-metabolite one-sided reaction such as
-    %"0.5 A + 0.5 B =>" by coincidence, when its magnitudes summed to
-    %exactly 1 --- and even then mishandled it below, since find() on such
-    %a reaction returns more than one index and only the first metabolite
-    %was ever used to build the boundary metabolite, silently dropping the
-    %rest.
+    %products. The coefficient magnitudes do not matter, so a scaled sink
+    %such as "2 A =>" and a multi-metabolite one-sided reaction such as
+    %"0.5 A + 0.5 B =>" are both boundary reactions.
     if any(col>0) && any(col<0)
         continue
     end
@@ -53,7 +47,7 @@
         continue
     end
     %A single metabolite names and inherits properties from that
-    %metabolite, as before. Several metabolites on the one populated side
+    %metabolite. Several metabolites on the one populated side
     %have no single metabolite to name the boundary after or copy
     %properties from, so it is named after the reaction instead and left
     %without those properties.

src/manipulation/setExchangeBounds.m

@@ -200,10 +200,8 @@
 
 % check for any metabolites that are exchanged in more than one reaction.
 % metInd holds one entry per (metabolite,reaction) pair, so a metabolite
-% involved in more than one exchange reaction is a repeated *value* in
-% metInd, not a position past numel(unique(metInd)) -- comparing those two
-% index spaces directly (as a previous version of this check did) compares
-% positions in the unique list against indices into metInd itself.
+% involved in more than one exchange reaction is a value that occurs more
+% than once in metInd. The occurrences are counted per unique value.
 [uniqueMetInd,~,ic] = unique(metInd);
 counts = accumarray(ic,1);
 multiMetInd = exchMetInd(uniqueMetInd(counts>1));

src/manipulation/simplifyModel.m

@@ -90,9 +90,10 @@
 end
 
 if deleteDuplicates==true
-    %Delete all but the last occurrence of duplicate reactions. The
-    %reactions must have the same bounds, reversibility, and objective
-    %coefficient to be regarded as duplicate
+    %Merge duplicate reactions into the first occurrence and delete the
+    %others. Reactions with the same stoichiometry and reversibility are
+    %duplicates; the merged reaction keeps the widest bounds and the largest
+    %objective coefficient (see contractModel)
     [reducedModel, rxnsToDelete, ~]=contractModel(reducedModel);
     deletedReactions=[deletedReactions; rxnsToDelete];
 end

src/queries/getIndexes.m

@@ -25,7 +25,8 @@
 % --------------------
 % returnLogical : logical
 %     sets whether to return a logical array or an array with the indexes
-%     (default false).
+%     (default false). It cannot be true for type 'metnames', which raises
+%     the error RAVEN:badInput.
 %
 % Returns
 % -------
@@ -61,6 +62,12 @@
     case 'genes'
         searchIn=model.genes;
     case 'metnames'
+        if returnLogical==true
+            error('RAVEN:badInput', ['The combination of type ''metnames'' and ' ...
+                'returnLogical true is not supported, because each name can ' ...
+                'match several metabolites. Use type ''mets'' or ''metcomps'' ' ...
+                'for a logical vector.']);
+        end
         searchIn=model.metNames;
     case 'metcomps'
         % If provided as metaboliteName[comp], then index search is quite

src/queries/printFluxes.m

@@ -46,8 +46,8 @@ function printFluxes(model, fluxes, varargin)
 %
 % Examples
 % --------
-%     printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile, ...
-%                 outputString, metaboliteList);
+%     printFluxes(model, fluxes, 'onlyExchange', false, 'cutOffFlux', 10^-6, ...
+%                 'outputFile', 'fluxes.txt', 'metaboliteList', {'glucose'});
 
 p=parseRAVENargs(varargin, {'onlyExchange',true; 'cutOffFlux',10^-8; 'outputFile',[]; 'outputString',[]; 'metaboliteList',{}});
 onlyExchange=p.onlyExchange;

src/reconstruction/kegg/getModelFromKEGG.m

@@ -24,7 +24,9 @@
 % Returns
 % -------
 % model : struct
-%     a model structure with all reactions and the metabolites used in them.
+%     a model structure with the reactions that the keep* options retain (by
+%     default all except those labelled as "general reaction") and the
+%     metabolites used in them.
 % KOModel : struct
 %     a model structure representing KEGG Orthology ids and their associated
 %     genes. KO ids are stored as reactions.

tests/function_tests/tAnalysis.m

@@ -122,6 +122,24 @@ function findGeneDeletionsComputesGrRatio(testCase)
             testCase.verifyGreaterThan(max(grRatioMuts), 0);
         end
 
+        function findGeneDeletionsDoubleDeletionSizes(testCase)
+            % Three genes with solvable single deletions give three gene pairs.
+            m = struct();
+            m.id = 'toy'; m.name = 'toy';
+            m.rxns = {'EX_a'; 'R1'; 'R2'; 'R3'; 'BIO'; 'EX_p'}; m.rxnNames = m.rxns;
+            m.mets = {'a'; 'b'; 'p'}; m.metNames = m.mets; m.metComps = ones(3,1);
+            m.comps = {'c'}; m.compNames = {'c'};
+            m.S = sparse([1 -1 -1 -1 0 0; 0 1 1 0 -1 0; 0 0 0 1 0 -1]);
+            m.lb = zeros(6,1); m.ub = [10; 1000; 1000; 1000; 1000; 1000];
+            m.rev = zeros(6,1); m.c = [0; 0; 0; 0; 1; 0]; m.b = zeros(3,1);
+            m.genes = {'g1'; 'g2'; 'g3'}; m.grRules = {''; 'g1'; 'g2'; 'g3'; ''; ''};
+            m.rxnGeneMat = sparse([0 0 0; 1 0 0; 0 1 0; 0 0 1; 0 0 0; 0 0 0]);
+            evalc('[genes, fluxes, ~, ~, grRatioMuts] = findGeneDeletions(m, ''dgd'');');
+            testCase.verifySize(genes, [3 2]);
+            testCase.verifySize(fluxes, [numel(m.rxns) 3]);
+            testCase.verifySize(grRatioMuts, [1 3]);
+        end
+
         function traceFluxPathReturnsCells(testCase)
             testCase.assumeSolver('solveLP');
             sol = solveLP(testCase.model);

tests/function_tests/tGapfilling.m

@@ -26,6 +26,28 @@ function checkProductionReturnsIndices(testCase)
             testCase.verifyClass(notProduced, 'double');
         end
 
+        function checkProductionReturnsExcretionReactionsOnly(testCase)
+            % Neither a nor c has a source, so b and d, and a and c, cannot be
+            % produced. The returned model holds an excretion reaction for
+            % each metabolite, with or without checkNeededForProduction.
+            m = struct();
+            m.id = 'toy'; m.name = 'toy';
+            m.rxns = {'R1'; 'R2'}; m.rxnNames = m.rxns;
+            m.mets = {'a'; 'b'; 'c'; 'd'}; m.metNames = m.mets; m.metComps = ones(4,1);
+            m.comps = {'c'}; m.compNames = {'c'};
+            m.S = sparse([-1 0; 1 0; 0 -1; 0 1]);
+            m.lb = [0; 0]; m.ub = [1000; 1000]; m.rev = [0; 0]; m.c = [0; 0];
+            m.b = zeros(4,1);
+            m.genes = {}; m.grRules = {''; ''}; m.rxnGeneMat = sparse(2,0);
+
+            [~, ~, ~, ~, plain] = checkProduction(m, 'printDetails', false);
+            [notProduced, ~, ~, ~, checked] = checkProduction(m, ...
+                'checkNeededForProduction', true, 'printDetails', false);
+            testCase.verifyEqual(notProduced, (1:4)');
+            testCase.verifyEqual(full(checked.S), full(plain.S));
+            testCase.verifyTrue(all(full(sum(checked.S(:, 3:end), 1)) < 0));
+        end
+
         function findLeakMetaboliteProduceReturnsSolution(testCase)
             evalc('[sol, metabolite] = findLeakMetabolite(testCase.model, ''produce'');');
             testCase.verifyNotEmpty(sol);

tests/function_tests/tGenomeData.m

@@ -77,6 +77,25 @@ function renameModelGenesRenamesGenesAndGrRules(testCase)
             testCase.verifyEqual(size(m2.rxnGeneMat,2), numel(m2.genes));
         end
 
+        function renameModelGenesReadsTableFromFilePath(testCase)
+            model.id         = 'test';
+            model.rxns       = {'r1'; 'r2'};
+            model.genes      = {'LT1'; 'LT2'};
+            model.grRules    = {'LT1'; 'LT1 or LT2'};
+            model.rxnGeneMat = sparse([1 0; 1 1]);
+
+            T = table({'LT1'; 'LT2'}, {'AAA'; 'BBB'}, ...
+                'VariableNames', {'locus_tag', 'gene_name'});
+            tsvFile = [tempname '.tsv'];
+            writetable(T, tsvFile, 'FileType', 'text', 'Delimiter', '\t');
+            testCase.addTeardown(@() delete(tsvFile));
+            output = evalc('m2 = renameModelGenes(model, tsvFile, ''locus_tag'', ''gene_name'');');
+
+            testCase.verifySubstring(output, tsvFile);
+            testCase.verifyEqual(sort(m2.genes), {'AAA'; 'BBB'});
+            testCase.verifyEqual(m2.grRules, {'AAA'; 'AAA or BBB'});
+        end
+
         function downloadGenomeDataDownloads(testCase)
             % Network test: download a small genome from NCBI Datasets.
             % Skipped when there is no internet connection.

tests/function_tests/tQueries.m

@@ -322,6 +322,19 @@ function getIndexesLogical(testCase)
             testCase.verifyEqual(find(lg), 1);
         end
 
+        function getIndexesMetnamesReturnLogicalErrors(testCase)
+            name = testCase.model.metNames(1);
+            testCase.verifyError( ...
+                @() getIndexes(testCase.model, name, 'metnames', 'returnLogical', true), ...
+                'RAVEN:badInput');
+            try
+                getIndexes(testCase.model, name, 'metnames', 'returnLogical', true);
+            catch e
+                testCase.verifySubstring(e.message, 'metnames');
+                testCase.verifySubstring(e.message, 'returnLogical');
+            end
+        end
+
         function getIndexesLogicalMaskAllTrueReturnsIndices(testCase)
             % A logical all-true mask of length n must return numeric [1..n],
             % not be passed through as a logical array: this requires an