Add various helpers: parseYAML, applyCondition, biomass, SBO, deltaG, findDuplicateRxns, curateModelFromTables

core/applyCondition.m

@@ -0,0 +1,130 @@
+function model = applyCondition(model, condition)
+% applyCondition
+%   Apply a deterministic "condition" to a model: a prelude that resets
+%   exchange bounds, optional metabolite removals + automatic charge
+%   rebalancing of a pseudoreaction, optional biomass-stoichiometry
+%   delta, and a per-reaction bounds diff. The schema is intentionally
+%   narrow so a condition can be reviewed as data.
+%
+%   Yeast-GEM was the first consumer; the same schema works for any
+%   GEM that keeps its condition presets as data rather than as code.
+%   Project-specific extensions (e.g. yeast-GEM's amino_acid_ratio
+%   step that rewrites a protein pseudoreaction's stoichiometry from a
+%   side-car TSV) are handled by the *caller* before / after this
+%   function — kept upstream-narrow on purpose.
+%
+%   Inputs:
+%       model       RAVEN model struct.
+%       condition   Either a path to a YAML condition file or a struct
+%                   already produced by parseYAML. The expected schema
+%                   (all keys optional):
+%
+%                       prelude:
+%                         reset_exchanges: out      % truthy -> reset all
+%
+%                       cofactor_pseudoreaction:
+%                         rxn_id: r_4598
+%                         remove_mets:
+%                           - { met: s_3714 }
+%                         charge_balance_met: s_0794
+%
+%                       biomass_stoichiometry_delta:
+%                         rxn_id: r_4041
+%                         add:
+%                           - { met: s_0689, coef:  0.08 }
+%                           - { met: s_0687, coef: -0.08 }
+%                           - { met: s_0794, coef: -0.16 }
+%
+%                       bounds:
+%                         - { rxn: r_1654, lb: -1000 }
+%                         - { rxn: r_1992, lb: 0 }
+%                         - { rxn: r_1663, lb: 0, ub: 0 }
+%
+%                       expected_uptake_count: 15
+%
+%   Output:
+%       model       Modified model.
+%
+% Usage: model = applyCondition(model, 'data/conditions/anaerobic.yml')
+%        model = applyCondition(model, parseYAML('data/conditions/anaerobic.yml'))
+
+if ischar(condition) || isstring(condition)
+    cond = parseYAML(char(condition));
+elseif isstruct(condition)
+    cond = condition;
+else
+    error('applyCondition:invalidCondition', ...
+        'condition must be a YAML file path or a struct.');
+end
+
+% --- Step 1: prelude ---------------------------------------------------
+if isfield(cond, 'prelude') && isfield(cond.prelude, 'reset_exchanges')
+    [~, exchangeRxns] = getExchangeRxns(model, cond.prelude.reset_exchanges);
+    model.lb(exchangeRxns) = 0;
+    model.ub(exchangeRxns) = 1000;
+end
+
+% --- Step 2: cofactor pseudoreaction edits ----------------------------
+if isfield(cond, 'cofactor_pseudoreaction')
+    cp = cond.cofactor_pseudoreaction;
+    cofacIdx = getIndexes(model, cp.rxn_id, 'rxns');
+    if isfield(cp, 'remove_mets')
+        for i = 1:numel(cp.remove_mets)
+            metIdx = getIndexes(model, cp.remove_mets{i}.met, 'mets');
+            model.S(metIdx, cofacIdx) = 0;
+        end
+    end
+    if isfield(cp, 'charge_balance_met')
+        balanceIdx = find(strcmp(model.mets, cp.charge_balance_met));
+        model.S(balanceIdx, cofacIdx) = 0;
+        model.S(balanceIdx, cofacIdx) = ...
+            -sum(model.S(:, cofacIdx) .* model.metCharges, 'omitnan');
+    end
+end
+
+% --- Step 3: biomass stoichiometry delta ------------------------------
+if isfield(cond, 'biomass_stoichiometry_delta')
+    delta = cond.biomass_stoichiometry_delta;
+    bioIdx = getIndexes(model, delta.rxn_id, 'rxns');
+    if isfield(delta, 'add')
+        for i = 1:numel(delta.add)
+            entry = delta.add{i};
+            metIdx = getIndexes(model, entry.met, 'mets');
+            model.S(metIdx, bioIdx) = full(model.S(metIdx, bioIdx)) + entry.coef;
+        end
+    end
+end
+
+% --- Step 4: bounds ---------------------------------------------------
+nUptake = 0;
+if isfield(cond, 'bounds')
+    for i = 1:numel(cond.bounds)
+        b = cond.bounds{i};
+        rxnIdx = find(strcmp(model.rxns, b.rxn));
+        if isempty(rxnIdx)
+            warning('applyCondition:missingRxn', ...
+                'Reaction %s not found in model; skipping.', b.rxn);
+            continue;
+        end
+        if isfield(b, 'lb')
+            model.lb(rxnIdx) = b.lb;
+            if b.lb == -1000
+                nUptake = nUptake + 1;
+            end
+        end
+        if isfield(b, 'ub')
+            model.ub(rxnIdx) = b.ub;
+        end
+    end
+end
+
+% --- Step 5: uptake sanity check --------------------------------------
+if isfield(cond, 'expected_uptake_count')
+    if nUptake ~= cond.expected_uptake_count
+        warning('applyCondition:uptakeCountMismatch', ...
+            'Expected %d uptake reactions, applied %d. Some may be missing from the model.', ...
+            cond.expected_uptake_count, nUptake);
+    end
+end
+
+end

core/assignSBOterms.m

@@ -0,0 +1,159 @@
+function model = assignSBOterms(model, opts)
+% assignSBOterms
+%   Assign SBO terms to metabolites and reactions following a generic
+%   rule set. Mirrors raven_python.annotation.add_sbo_terms;
+%   organism-agnostic, parameterised entirely by `opts`. The
+%   yeast-GEM port of this function is the legacy addSBOterms.m,
+%   which becomes a thin shim here.
+%
+%   Rules
+%   -----
+%   Metabolites:
+%       SBO:0000649 (Biomass) when met.name is in opts.biomassMetNames,
+%       or ends with any of opts.biomassMetSuffixes. Otherwise
+%       SBO:0000247 (Simple chemical).
+%
+%   Reactions (default → override → pseudoreaction override):
+%       SBO:0000176 (Metabolic reaction) default.
+%       Single-reactant reactions become:
+%           SBO:0000627 (exchange) if the lone metabolite is
+%               extracellular (compartment 'e' or compartment name
+%               containing 'extracellular'),
+%           SBO:0000632 (sink) if coef < 0,
+%           SBO:0000628 (demand) otherwise.
+%       Transport reactions (detected by opts.transportDetector or
+%           the default heuristic: same metName in ≥ 2 compartments
+%           in a single reaction) → SBO:0000655.
+%       Reactions whose name matches opts.biomassRxnName → SBO:0000629.
+%       Reactions whose name matches opts.ngamRxnName → SBO:0000630.
+%       Reactions whose name contains any of
+%           opts.pseudoreactionSubstrings → SBO:0000395.
+%
+%   "fill" semantic — SBO is written via editMiriam(..., 'fill') so
+%   pre-existing SBO annotations are preserved.
+%
+%   Inputs:
+%       model   RAVEN model struct.
+%       opts    (opt) struct with any of the following fields. Missing
+%               fields take the defaults shown:
+%                 biomassMetNames        {'biomass','DNA','RNA','protein',
+%                                         'carbohydrate','lipid','cofactor','ion'}
+%                 biomassMetSuffixes     {' backbone',' chain'}
+%                 biomassRxnName         'biomass pseudoreaction'
+%                 ngamRxnName            'non-growth associated maintenance reaction'
+%                 pseudoreactionSubstrings {'pseudoreaction','SLIME rxn'}
+%                 onlyLastReactionForPseudo  false. yeast-GEM bug-compat
+%                                            flag — replicates the
+%                                            legacy `for i=numel(...)`
+%                                            typo (pseudoreaction
+%                                            SBOs applied only to the
+%                                            last reaction). Off by
+%                                            default; turn ON for
+%                                            byte-equivalent yeast-GEM
+%                                            output.
+%
+%   Output:
+%       model   Modified model.
+%
+% Usage: model = assignSBOterms(model)
+%        model = assignSBOterms(model, struct('onlyLastReactionForPseudo', true))
+
+if nargin < 2 || isempty(opts)
+    opts = struct();
+end
+opts = applyDefaults(opts);
+
+% Metabolite SBO ------------------------------------------------------
+metsSBO = cell(size(model.mets));
+for i = 1:length(model.mets)
+    metName = model.metNames{i};
+    if any(strcmp(opts.biomassMetNames, metName)) || endsWithAny(metName, opts.biomassMetSuffixes)
+        metsSBO{i} = 'SBO:0000649';
+    else
+        metsSBO{i} = 'SBO:0000247';
+    end
+end
+
+% Reaction SBO --------------------------------------------------------
+rxnSBO = cell(size(model.rxns));
+rxnSBO(:) = {'SBO:0000176'};
+
+% Single-reactant reactions
+reactantNumber = sum(model.S ~= 0, 1);
+singleRxns = find(reactantNumber == 1);
+for k = 1:numel(singleRxns)
+    idx = singleRxns(k);
+    metRow = find(model.S(:, idx));
+    compName = model.compNames{model.metComps(metRow)};
+    compShort = model.comps{model.metComps(metRow)};
+    if strcmp(compShort, 'e') || strcmp(compName, 'extracellular')
+        rxnSBO{idx} = 'SBO:0000627';
+    elseif sum(model.S(:, idx)) < 0
+        rxnSBO{idx} = 'SBO:0000632';
+    else
+        rxnSBO{idx} = 'SBO:0000628';
+    end
+end
+
+% Transport reactions
+if isfield(opts, 'transportRxnIdxs') && ~isempty(opts.transportRxnIdxs)
+    transportIdxs = opts.transportRxnIdxs;
+else
+    transportIdxs = getTransportRxns(model);
+end
+rxnSBO(transportIdxs) = {'SBO:0000655'};
+
+% Pseudoreaction overrides
+if opts.onlyLastReactionForPseudo
+    pseudoTargets = numel(model.rxns);
+else
+    pseudoTargets = 1:numel(model.rxns);
+end
+for ii = pseudoTargets
+    name = model.rxnNames{ii};
+    if strcmp(name, opts.biomassRxnName)
+        rxnSBO{ii} = 'SBO:0000629';
+    elseif strcmp(name, opts.ngamRxnName)
+        rxnSBO{ii} = 'SBO:0000630';
+    else
+        for k = 1:numel(opts.pseudoreactionSubstrings)
+            if contains(name, opts.pseudoreactionSubstrings{k})
+                rxnSBO{ii} = 'SBO:0000395';
+                break;
+            end
+        end
+    end
+end
+
+model = editMiriam(model, 'met', 'all', 'sbo', metsSBO, 'fill');
+model = editMiriam(model, 'rxn', 'all', 'sbo', rxnSBO, 'fill');
+end
+
+
+function opts = applyDefaults(opts)
+defaults = struct( ...
+    'biomassMetNames', {{'biomass','DNA','RNA','protein','carbohydrate','lipid','cofactor','ion'}}, ...
+    'biomassMetSuffixes', {{' backbone',' chain'}}, ...
+    'biomassRxnName', 'biomass pseudoreaction', ...
+    'ngamRxnName', 'non-growth associated maintenance reaction', ...
+    'pseudoreactionSubstrings', {{'pseudoreaction','SLIME rxn'}}, ...
+    'onlyLastReactionForPseudo', false);
+fields = fieldnames(defaults);
+for k = 1:numel(fields)
+    f = fields{k};
+    if ~isfield(opts, f) || isempty(opts.(f))
+        opts.(f) = defaults.(f);
+    end
+end
+end
+
+
+function tf = endsWithAny(s, suffixes)
+tf = false;
+for i = 1:numel(suffixes)
+    if endsWith(s, suffixes{i})
+        tf = true;
+        return;
+    end
+end
+end