Skip to content

fix: align raven-toolbox with RAVEN (parity review, round 3) - #186

Merged
edkerk merged 40 commits into
developfrom
parity/round-3
Oct 30, 2026
Merged

edkerk merged 40 commits into
developfrom
parity/round-3

Conversation

@edkerk

@edkerk edkerk commented Oct 30, 2026 •

Copy link
Copy Markdown
Collaborator

Third round of the behaviour parity review of raven-toolbox against RAVEN. The matching RAVEN changes are in CompMetEng/RAVEN#15 and the scenarios and ledger rows in CompMetEng/raven-gecko-parity#143; the previous round is #182.

Defaults and behaviour

Function Before After
manipulation.remove_genes default blocked_reactions "remove" "constrain", as removeGenes: a reaction whose rule loses its last enzyme stays, with the rule dropped and bounds (0, 0)
manipulation.remove_genes genes the requested genes were removed when at least one requested gene is in the model, every gene that no remaining reaction named before the removal is removed as well; a gene kept by a constrained reaction stays; when none is in the model, nothing changes
manipulation.merge_compartments surviving id first member with its suffix rewritten to the merged compartment (ACCOA_s) the copy already in the merged compartment, otherwise the id that sorts first (ACCOA_c for ACCOA_c and ACCOA_m), as mergeCompartments
manipulation.merge_compartments annotation dropped kept from the surviving copy
manipulation.merge_compartments default drop_single_metabolite_reactions True False, as deleteRxnsWithOneMet
default metabolite key of localization.assign_compartments, apply_assignment, relocate_reactions, curation_priority the id without the suffix of the metabolite's own compartment notes["merged_base_id"] when merge_compartments stored it, otherwise as before, so a merged ACCOA_c in s and a universal ACCOA_m share the key ACCOA
two metabolites with one key in one compartment (same four functions) one of them stood for both ValueError
default threads of reconstruction.homology.run_blast, run_diamond, run_supporth and reconstruction.kegg.run_hmmsearch, get_kegg_model_from_sequences, build_ko_hmm, build_hmm_library the machine's CPU count minus one every CPU the process may run on, as getBlast, getDiamond, runSuppOrth and getKEGGModelForOrganism
gapfilling.analyse_topology default seeds one-metabolite reactions with lower_bound < 0 metabolites of every reaction without products and with lower_bound < 0, or without substrates and with upper_bound > 0; UserWarning when there is none
gapfilling.analyse_topology default targets every metabolite of the model substrates and products of the first objective reaction; ValueError without an objective
init.merge_linear group numbers merge order order of the first reaction of each group in the original model, as mergeLinear
io.read_yaml_model notes written as a notes mapping nested under notes["note"] back in notes

Gene rules

New utils.parse_gr_rule reads a rule in the grammar of RAVEN's parseGrRule (and/or in any case, &, &&, |, ||, operators as whole words, whitespace as MATLAB's isspace), so G1&G2 is one gene id; a malformed rule raises ValueError with the message of RAVEN:badGrRule. manipulation.gpr_to_dnf and utils.is_dnf take a rule string and read it with this parser. utils.gpr.set_gene_rule sets a rule from a GPR object or a string read with parse_gr_rule, and io.read_yaml_model, curation.batch_curate (and batch_curate_from_tsv) and the functions that copy or build a rule (manipulation.merge_compartments, copy_to_compartment, convert_to_irreversible, merge_models, add_reactions_from_model, expand_model, change_gene_reaction_rules, add_reactions_from_equations, localization.apply_assignment, init.fill_tasks, init.remove_low_score_genes, reconstruction.homology.get_model_from_homology, reconstruction.kegg.assemble_model_from_ko_genes, io.read_matlab_model) set rules through it, so G1&G2 or G3 has the genes G1&G2 and G3, as in RAVEN; a rule that parse_gr_rule refuses is read by cobra. cobra's JSON reader still splits G1&G2.

merge_compartments stores the base id of each merged metabolite in notes["merged_base_id"] (ACCOA for the kept ACCOA_c, the own id when two would coincide). The note survives raven-toolbox YAML, copy(), cobra JSON and SBML; RAVEN's readYAMLmodel shows it in metNotes and no RAVEN function reads it, and a round trip through RAVEN YAML or a .mat file loses it, after which the localization functions key a merged ACCOA_c in s as ACCOA_c without an error.

ftINIT tie resolution

With resolve_ties, init.ftinit and init.run_ftinit take each step's fluxes from the least-total-flux solution for the tie-broken on set (phase 2c, gap 1e-6, full time_limit) and, among those, from the one with the least id-rank-weighted backward flux (phase 2d), so the directions the next step forces no longer depend on the solver. The id-rank phases rank a merged reaction by the smallest original reaction id in its group (new run_ftinit argument rank_keys), and a merged reaction prefers its forward direction as written. The last step of ftinit skips phases 2c and 2d, as no step reads its fluxes, and reports the fluxes of the id-rank phase; the reaction set does not change. A phase that finds no solution raises a UserWarning naming it (phase2c-minflux, phase2d-direction) and the fluxes of the phase before it are kept; a phase that stops at the time limit is adopted with a UserWarning.

With CompMetEng/RAVEN#15, the 88 RAVEN scenarios of CompMetEng/raven-gecko-parity#143 give 86 MATCH and 2 DIFFER, both documented in its ledger: in fseof_toy, solveLP relaxes EX_S to 0.999999 of its maximum, and in edge_model_toy, writeYAMLmodel refuses two metabolites with one name in one compartment where raven-toolbox warns and writes the file. On Human-GEM with HCT116 (series 1+1), both tools find 204 task-essential reactions and give identical final reaction sets (7,752 reactions); 28 forced directions of step 2 differ, as the ledger records. For Human-GEM's gene-essentiality baseline, the template-bounds output already on develop (c3a2130) means the complex III and IV genes no longer block growth in the growth-ratio score (139 growth-lethal flips summed over five lines against about 18 from run-to-run noise), matching RAVEN, because Human-GEM's own succinate dehydrogenase direction is reversible; viability and task MCCs are unchanged within noise.

Not checked: GECKO scenarios and geckopy, CI on Python 3.11 to 3.15, RAVEN test classes outside tUtils, tManipulation, tINIT, tGapfilling, tAnalysis and tReconstruction, gene essentiality on Human-GEM beyond one run per variant.

edkerk added 30 commits October 28, 2026 15:11
… solution

After the id-rank phase, a further phase holds the on/off state and minimises the summed absolute flux, so the direction and magnitude at which the next step forces a reaction do not depend on the solver's vertex.
… ranks a merged reaction by its smallest original id
… or the one that sorts first

The merged metabolite took the first member's id with its compartment
suffix rewritten (ACCOA_c became ACCOA_s). It now keeps the id of the copy
already in the merged compartment, if any, otherwise the id that sorts
first, as RAVEN's mergeCompartments does.
Group numbers followed the merge order. They are now 1, 2, ... in the
order of the first reaction of each group in the original model, as
RAVEN's mergeLinear numbers them.
…mar of RAVEN's parseGrRule, and raise ValueError for a malformed rule
…_metabolite_reactions defaults to False

The merged metabolite kept the id, name, formula, charge and notes of the
chosen copy and dropped its annotation. It now keeps the annotation too, as
RAVEN's mergeCompartments keeps every field of the kept row.
drop_single_metabolite_reactions defaults to False, the default of RAVEN's
deleteRxnsWithOneMet.
…e universal copy

merge_compartments keeps the id of one copy, so a merged ACCOA_c sits in s and
its own-compartment suffix rule left the key ACCOA_c. A universal ACCOA_m has
the key ACCOA, so apply_assignment added a second copy of the metabolite in m.
assign_compartments, apply_assignment, relocate_reactions and
curation_priority now also strip a suffix that names a compartment of the
model, the universal model, the scores or the placements.
The tokeniser split words at Python's Unicode whitespace, so a no-break space
(U+00A0), U+0085, U+2007 and U+202F separated words, while RAVEN's
parseGrRule, which uses isspace, keeps them in the word and refuses
G1<U+00A0>and G2. The whitespace class is now the list of code points for
which isspace(char(0:65535)) is true in MATLAB R2024b.
The rule was checked with parse_gr_rule, for which G1&G2 is one gene id, and
then stored as a string that cobra reads as G1 and G2. The reaction now gets
the parsed rule, so its genes are the ones that were checked, as in RAVEN's
curateModelFromTables.
…y CPU the process may use

The default was the machine's CPU count minus one. It is now the number of
CPUs the process may run on, as in run_blast and run_diamond and in RAVEN's
getKEGGModelForOrganism, which passes feature('numcores')'s assigned count
to hmmsearch.
…g go back into notes

write_yaml_model writes the notes without a RAVEN field as a notes mapping,
which read_yaml_model nested under notes["note"]. They are now spread back
into notes, with the free-text note under "note".
…ments stores

The default key guessed a compartment suffix from the compartments each
function saw, so assign_compartments and apply_assignment could key one
metabolite differently (I and I_g), and distinct species such as X_c and
X_m in s, or vitamin_a and vitamin_c, got one key. merge_compartments now
stores the base id of each merged metabolite in notes["merged_base_id"]
(the own id when two would coincide), and the default key is that stored id,
otherwise the id without the suffix of the metabolite's own compartment.
assign_compartments, apply_assignment, relocate_reactions and
curation_priority raise ValueError when two metabolites in one compartment
have the same key.
cobra's model_from_dict read the rule G1&G2 or G3, written by write_yaml_model
after batch_curate, as the genes G1, G2 and G3. The reader now gives each
reaction the rule as parse_gr_rule reads it, so G1&G2 stays one gene; a rule
that parse_gr_rule refuses is read by cobra as before. The batch_curate help
states that cobra's JSON reader still splits such an id.
…ry CPU the process may use

As run_blast and run_hmmsearch, and as RAVEN's runSuppOrth, which passes the
cores assigned to MATLAB (was: the machine's CPU count minus one).
The test patched a helper that only exists once the default is in place, so
without it the test stopped with an AttributeError. It now fixes the CPUs
the process may use and checks the --cpu value, also for
get_kegg_model_from_sequences without threads.
utils.gpr.set_gene_rule copies a GPR object, or reads a rule string with
parse_gr_rule, so G1&G2 stays one gene id. merge_compartments,
copy_to_compartment, apply_assignment, convert_to_irreversible,
merge_models, add_reactions_from_model, expand_model,
change_gene_reaction_rules, add_reactions_from_equations, fill_tasks,
remove_low_score_genes, the homology and KEGG model builders and
read_matlab_model set rules through it instead of assigning the string
to gene_reaction_rule, which cobra read with its own grammar.
… written

_min_backward_flux weights the backward flux by the rank key only and a turned-on reaction with zero flux is forced forward in the next step. run_ftinit takes the private keyword _pin_fluxes, which ftinit sets false in its last step. The tests check that run_ftinit skips phases 2c and 2d with _pin_fluxes false and that ftinit sets it false in the last step only.
@edkerk
edkerk merged commit 94f3940 into develop Oct 30, 2026
10 checks passed
@edkerk
edkerk deleted the parity/round-3 branch October 30, 2026 11:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant