Fixes from the RAVEN 3 Guide review

src/analysis/findGeneDeletions.m

@@ -36,8 +36,10 @@
 %     with details about each gene in originalGenes and why or why not it
 %     was deleted:
 %
-%     - 1 : was deleted
-%     - 2 : proved lethal in sgd (single gene deletion)
+%     - 1 : was deleted, and the problem was solved (a lethal deletion
+%       with zero growth also gets 1)
+%     - 2 : the single-deletion problem had no optimal solution
+%       (infeasible or failed); fluxes and grRatioMuts are 0
 %     - 3 : redundant, no longer used
 %     - 4 : involved in dead-end reaction
 % grRatioMuts : double

src/analysis/getMinimalMedium.m

@@ -50,6 +50,16 @@
 %
 % A y_i = 0 forces reaction i to carry non-negative flux (no net uptake);
 % y_i = 1 allows it to carry flux as low as lb_i.
+%
+% The minimal number of uptakes is guaranteed, but the medium itself need
+% not be unique: when several media of the same size support the growth
+% rate, the solver returns one of them, and another solver or solver
+% version can return a different one.
+%
+% The MILP is solved with the solver set by setRavenSolver. RAVEN's bundled
+% GLPK solves it for small models (with a RAVEN:glpkMILP warning); for
+% larger models Gurobi is much faster and more reliable
+% (setRavenSolver('gurobi')).
 
 p = parseRAVENargs(varargin, {'minGrowth',[]; 'verbose',true});
 minGrowth = p.minGrowth;

src/analysis/haveFlux.m

@@ -1,9 +1,10 @@
 function I=haveFlux(model,varargin)
 % haveFlux  Check which reactions can carry a flux.
 %
-% Checks which reactions can carry a (positive or negative) flux. Is used as
-% a faster version of getAllowedBounds if it is only interesting whether the
-% reactions can carry a flux or not.
+% Checks which reactions can carry a (positive or negative) flux under the
+% model's current bounds; other media or bounds are not explored. Is used
+% as a faster version of getAllowedBounds if it is only interesting whether
+% the reactions can carry a flux or not.
 %
 % Parameters
 % ----------
@@ -28,7 +29,7 @@
 % -------
 % I : logical
 %     logical array with true if the corresponding reaction can carry a
-%     flux.
+%     flux under the tested bounds.
 %
 % Examples
 % --------
@@ -38,6 +39,12 @@
 % -----
 % If a model has +/- Inf bounds then those are replaced with an arbitary
 % large value of +/- 10000 prior to solving.
+%
+% Dead-end reactions are removed (simplifyModel) before solving. A dead-end
+% reaction with a forced flux (lb > 0 or ub < 0) makes the model infeasible,
+% but it is removed like any other dead end, so the remaining reactions are
+% tested as if it were absent. A forced flux elsewhere in an infeasible
+% model gives false for all reactions.
 
 p=parseRAVENargs(varargin, {'cutOff',10^-6; 'rxns',[]; 'seed',[]});
 if ~isempty(p.seed)
@@ -66,6 +73,13 @@
 %First remove all dead-end reactions
 smallModel=simplifyModel(model,false,false,true,true);
 
+%Nothing is left to solve when every reaction is a dead end, and reactions
+%without metabolites are limited by their bounds only
+if isempty(smallModel.rxns) || isempty(smallModel.mets)
+    I=ismember(rxns,smallModel.rxns(smallModel.ub>cutOff | smallModel.lb<-cutOff));
+    return
+end
+
 %Get the indexes of the reactions in the connected model
 indexes=getIndexes(smallModel,intersect(smallModel.rxns,rxns),'rxns');
 J=false(numel(indexes),1);

src/analysis/randomSampling.m

@@ -81,9 +81,12 @@
 %
 % Notes
 % -----
-% The solutions are generated by maximizing (with random weights) for a
-% random set of nObjectives reactions (default 2). For reversible reactions
-% it randomly chooses between maximizing and minimizing.
+% With 'randomObjective', the solutions are generated by maximizing (with
+% random weights) for a random set of nObjectives reactions (default 2).
+% For reversible reactions it randomly chooses between maximizing and
+% minimizing. The default 'achr' and the 'chrr' methods use no random
+% objectives; they take hit-and-run steps through the interior of the
+% solution space.
 %
 % If the model is a GECKO v3+ ecModel, then usage_prot reactions are not
 % selected for sampling, instead focusing on sampling from the metabolic

src/analysis/reporterMetabolites.m

@@ -36,10 +36,10 @@
 %       supplied.
 %     - mets : a cell array of metabolite IDs for the metabolites for which
 %       a score could be calculated.
-%     - metZScores : Z-scores for differential expression around each
-%       metabolite in "mets".
-%     - metPValues : P-values for differential expression around each
-%       metabolite in "mets".
+%     - metZScores : background-corrected Z-scores for differential
+%       expression around each metabolite in "mets". The entries of all
+%       fields are sorted by this score, highest first.
+%     - metPValues : one-sided (upper-tail) P-values of "metZScores".
 %     - metNGenes : number of neighbouring genes for each metabolite in
 %       "mets".
 %     - meanZ : average Z-scores for the genes around each metabolite in
@@ -57,6 +57,10 @@
 % For details about the algorithm, see Patil KR, Nielsen J, Uncovering
 % transcriptional regulation of metabolism by using metabolic network
 % topology. Proc. Natl Acad. Sci. USA 2005;102:2685-2689.
+%
+% The background correction uses the closed-form mean and standard
+% deviation of the aggregate Z-score of n genes drawn at random from all
+% scored genes, so the results are deterministic.
 
 genes=convertCharArray(genes);
 p=parseRAVENargs(varargin, {'printResults',false; 'outputFile',[]; 'geneFoldChanges',[]});
@@ -136,11 +140,11 @@
 
 %Correct for background using the closed-form moments of the sampling
 %distribution. For n genes drawn with replacement from the scored-gene pool
-%(mean mu, std sigma), the aggregate Z = sum(z)/sqrt(n) has:
+%(mean mu, population std sigma), the aggregate Z = sum(z)/sqrt(n) has:
 %  E[Z] = sqrt(n)*mu  and  Std[Z] = sigma  (independent of n)
 %so the corrected score is (metZ - sqrt(n)*mu) / sigma.
 mu = mean(geneZScores);
-sigma = std(geneZScores);
+sigma = std(geneZScores,1);
 metZScores = (metZScores - sqrt(metNGenes) .* mu) ./ sigma;
 
 %Calculate the P-values (upper-tail of the standard normal). Uses the base

src/analysis/sampleChebyshevCenter.m

@@ -26,7 +26,7 @@
 norms = sqrt(sum(A.^2, 2));
 
 % Variables [x (n, free); r (1, >=0)]; minimise -r (maximises the ball radius).
-prob.A      = [A, norms];
+prob.A      = sparse([A, norms]);      % Gurobi and SCIP require a sparse A
 prob.a      = prob.A;
 prob.b      = b;
 prob.csense = repmat('L', 1, m);

src/annotation/assignSBOterms.m

@@ -6,7 +6,10 @@
 % organism-agnostic, parameterised entirely by opts.
 %
 % SBO is written via editMiriam(..., 'fill') so pre-existing SBO
-% annotations are preserved.
+% annotations are preserved. The terms are stored in model.metMiriams and
+% model.rxnMiriams under the namespace 'sbo' (there is no separate SBO
+% field); read them with extractMiriam(model.rxnMiriams, 'sbo').
+% exportModel writes them as the SBML sboTerm attribute.
 %
 % Parameters
 % ----------

src/biomass/scaleBiomassFraction.m

@@ -2,7 +2,8 @@
 % scaleBiomassFraction  Rescale a biomass component to a target value.
 %
 % Rescale a biomass component to a target g/gDW value, optionally
-% balancing a second component so the total biomass mass stays at 1 g/gDW.
+% rescaling a second component so that the biomass total becomes 1 g/gDW,
+% whatever the total was before.
 % The MATLAB counterpart of yeast-GEM's legacy scaleBioMass.
 %
 % Parameters
@@ -19,8 +20,10 @@
 % Name-Value Arguments
 % --------------------
 % balanceOut : char
-%     Second component name to adjust so the biomass total remains 1
-%     g/gDW. Empty / omit to skip balancing.
+%     Second component name to adjust so the biomass total becomes 1
+%     g/gDW. It must be large enough to absorb the change: an error is
+%     raised if its fraction would become zero or negative, or if it is
+%     0 to begin with. Empty / omit to skip balancing.
 %
 % Returns
 % -------
@@ -55,11 +58,24 @@
 if ~isempty(balanceOut)
     fractions = getBiomassFractions(model, biomassConfig);
     total = fractions.total;
+    if ~isfield(fractions, balanceOut)
+        error('scaleBiomassFraction:unknownComponent', ...
+            'biomassConfig has no component named %s', balanceOut);
+    end
     balanceCurrent = fractions.(balanceOut);
     if balanceCurrent == 0
-        return;
+        error('scaleBiomassFraction:zeroCurrent', ...
+            ['Cannot balance with %s: current fraction is 0 ' ...
+             '(pseudoreaction missing or empty).'], balanceOut);
+    end
+    balanceNew = balanceCurrent + (1 - total);
+    if balanceNew <= 0
+        error('scaleBiomassFraction:negativeBalance', ...
+            ['Cannot balance with %s: its fraction would become %g g/gDW. ' ...
+             'It is %g g/gDW, but the other components total %g g/gDW.'], ...
+            balanceOut, balanceNew, balanceCurrent, total - balanceCurrent);
     end
-    balanceFactor = (balanceCurrent + (1 - total)) / balanceCurrent;
+    balanceFactor = balanceNew / balanceCurrent;
     model = scaleBiomassPseudoreaction(model, biomassConfig, balanceOut, balanceFactor);
 end
 end