Reconstruction (Python)#

raven-toolbox objects in raven_toolbox.reconstruction.

Functions#

Function	Summary
`assemble_model_from_ko_genes`	Build a draft model from a `{ko: [gene, ...]}` assignment.
`assign_kos`	Assign genes to KOs from HMM hits, applying the cut-off and ratio filters.
`blast_from_table`	Load a precomputed homology hits table (CSV path or DataFrame, not Excel).
`build_hmm_library`	Build a domain (`"prokaryotes"`/`"eukaryotes"`) HMM library.
`build_kegg_tables`	Build the minimal relational tables from parsed records.
`build_ko_fastas`	Write one `<KO>.fa` per KO with its member genes' sequences.
`build_ko_hmm`	Cluster, align and train a profile HMM for one KO's multi-FASTA.
`build_reference_model`	Assemble the gene-free KEGG reference model from parsed records.
`download_kegg_dump`	Fetch and arrange a complete KEGG dump into `dest`.
`extract_kegg_dump`	Extract and arrange the downloaded archives into the flat dump layout.
`fetch_kegg_files`	Download the raw KEGG archives into `dest` (basenames). Returns the paths.
`flag_set`	Reaction ids whose `column` flag is truthy (handles bool or TSV strings).
`get_kegg_model_for_organism`	Reconstruct a draft model for a KEGG species from its KO annotations.
`get_kegg_model_for_organism_from_artefacts`	Load the published 3b.2 artefacts from `artefact_dir` and build the model.
`get_kegg_model_from_sequences`	Reconstruct a draft model for a proteome by HMM-searching the KO library.
`get_kegg_model_from_sequences_with_artefacts`	Load reference model + tables from `artefact_dir` and run the HMM query.
`get_model_from_homology`	Build a draft model for `model_for` by transferring reactions from templates.
`hits_from_supporth`	Build a homology hits table from a SuppOrth `suppOrth.tsv`.
`HomologyResult`	Result of `get_model_from_homology`.
`KeggCompound`	A metabolite parsed from the KEGG `compound` flat file.
`KeggKO`	A KEGG Orthology entry: its name and the organism genes assigned to it.
`KeggReaction`	A reaction parsed from the KEGG `reaction` flat file.
`make_ortholog_hits`	Build a bidirectional hits table from a predefined ortholog list.
`organism_domains`	Return `{organism_id: domain}` (the outermost category).
`organisms_in_domain`	Organism ids whose outermost category matches `domain` (case-insensitive).
`parse_hmmsearch_tblout`	Parse `hmmsearch --tblout` text into a `[ko, gene, evalue]` table.
`parse_kegg_compounds`	Parse `<kegg_dir>/compound` (+ optional `compound.inchi`) into records.
`parse_kegg_dump`	Parse a full KEGG dump into the reference model + tables and write them out.
`parse_kegg_kos`	Parse `<kegg_dir>/ko` into `KeggKO` records (name + organism genes).
`parse_kegg_reactions`	Parse `<kegg_dir>/reaction` into `KeggReaction` records.
`parse_taxonomy`	Return `{organism_id: [category, ...]}` from outermost to innermost.
`parse_taxonomy_records`	Return `[(organism_id, name, [category, ...]), ...]` in file order.
`phyl_dist`	Compute the KEGG phylogenetic distance matrix from the `taxonomy` file.
`PhylDist`	KEGG phylogenetic distance structure — RAVEN `getPhylDist` / `keggPhylDist`.
`read_kegg_table`	Read a KEGG table written by `write_kegg_tables` or
`run_blast`	Bidirectional BLAST+ between an organism and template organisms.
`run_diamond`	Bidirectional DIAMOND between an organism and template organisms.
`run_hmmsearch`	Search the profile `library` against proteome `fasta`; return tblout text.
`run_supporth`	Run SuppOrth on an organism and template proteomes and read its consensus.
`stream_organism_gene_ko`	Stream the `ko` file to a sorted, gzipped `organism_gene_ko.tsv.gz`.
`validate_hits`	Check a hits DataFrame has the required columns; return it unchanged.
`write_kegg_tables`	Write each table as a gzipped TSV (`<prefix><name>.tsv.gz`) into `out_dir`.

Reference#

assemble_model_from_ko_genes#

assemble_model_from_ko_genes(reference_model: cobra.Model, ko_reaction: pd.DataFrame, ko_to_genes: dict[str, list[str]], *, rxn_flags: pd.DataFrame | None=None, keep_spontaneous: bool=True, keep_undefined_stoich: bool=True, keep_incomplete: bool=True, keep_general: bool=False, model_id: str | None=None, model_name: str | None=None, note: str | None=None, prune_orthology: bool=False)

Build a draft model from a {ko: [gene, ...]} assignment.

Returns (model, gpr_map) where gpr_map is the kept reactions' gene lists, so callers can add gene annotations afterwards.

rxn_flags is a per-reaction quality-flag table (a reaction column plus any of spontaneous/undefined_stoich/incomplete/general). A reaction flagged undefined_stoich, incomplete or general is dropped unless its matching keep_* argument is True — this takes precedence over the reaction having genes. keep_spontaneous (default True) additionally keeps otherwise gene-less reactions flagged spontaneous.

prune_orthology restricts each gene-backed reaction's kegg.orthology annotation to the KOs that actually contributed a gene, matching RAVEN getKEGGModelForOrganism's HMM branch (the organism-annotation path keeps the full reference KO list, so this is off by default).

assign_kos#

assign_kos(hits: pd.DataFrame, *, cutoff: float=1e-30, min_score_ratio_ko: float=0.3, min_score_ratio_g: float=0.9)

Assign genes to KOs from HMM hits, applying the cut-off and ratio filters.

Ports RAVEN's three steps on the KO×gene E-value matrix:

keep hits with evalue <= cutoff;
min_score_ratio_ko — within a KO, drop genes whose log(evalue)/log(best_evalue_in_KO) < min_score_ratio_ko (prune weak members of a KO);
min_score_ratio_g — within a gene, drop KOs whose log(evalue)/log(best_evalue_for_gene) < min_score_ratio_g (stop a gene that clearly belongs to one KO leaking into weaker ones).

Smaller E-value = better; since all kept values are < 1 their logs are negative, so the best (smallest) hit gives ratio 1 and weaker hits give a smaller positive ratio.

Default calibration, cross-validated against the true KEGG gene→KO annotation of four organisms spanning the prok/euk libraries and the well-/lesser-studied axis (S. cerevisiae, Cyanidioschyzon merolae, E. coli, Mycoplasma genitalium): real annotations score overwhelmingly (median E ≈ 1e-100…1e-155) while spurious hits pile up at ≈1e-8, so the two are separated by ~20 orders of magnitude. RAVEN's 1e-50 sits inside the true tail and silently drops real but divergent hits — costing 16% gene→KO recall on the divergent minimal genome (M. genitalium) for no noise-rejection benefit (noise is far weaker). The default is therefore loosened to 1e-30 (recovers that tail; still ~22 orders above the noise floor), with the precision work moved to min_score_ratio_g = 0.9 — the effective precision lever (it resolves multi-KO genes). min_score_ratio_ko proved empirically inert across all four organisms (identical output at 0.0/0.3/0.5) and is kept only for RAVEN parity.

blast_from_table#

blast_from_table(source: str | Path | pd.DataFrame)

Load a precomputed homology hits table (CSV path or DataFrame, not Excel).

Must contain the HIT_COLUMNS columns.

build_hmm_library#

Build a domain ("prokaryotes"/"eukaryotes") HMM library.

Restricts genes to the domain's organisms (from taxonomy), builds a multi-FASTA and a profile HMM per KO under out_dir, and (if concatenate) concatenates them into out_dir/library.hmm — searched directly with hmmsearch (no hmmpress needed). Returns {"hmms": [...], "library": path | None}.

Heavy and binary-dependent — intended for the maintainer, run once per KEGG release. Resumable: skips KOs that already have an .hmm and reuses an already-built multi-FASTA (see build_ko_fastas), so re-running after a failure continues rather than restarting. force rebuilds everything.

progress shows tqdm bars for the multi-FASTA and per-KO HMM passes (the headline N of M counter for the long build); it is independent of verbose (per-KO INFO logging). With both on, per-KO log lines are routed through tqdm.write so they don't corrupt the bar.

build_kegg_tables#

build_kegg_tables(reactions: list[KeggReaction], kos: list[KeggKO])

Build the minimal relational tables from parsed records.

Returns a dict of DataFrames keyed by table name: ko_reaction, ko_names, organism_gene_ko, rxn_flags.

build_ko_fastas#

build_ko_fastas(organism_gene_ko: pd.DataFrame, genes_pep: str | Path, out_dir: str | Path, *, organisms: set[str] | None=None, progress: bool=False, force: bool=False)

Write one <KO>.fa per KO with its member genes' sequences.

Uses a stdlib offset index instead of a Java-hashtable byte scan. organisms restricts to a domain's organism codes (for the prok/euk split). Empty KOs are skipped (no file). progress shows a tqdm bar over the per-KO writing pass. Returns {ko: path} for the files written.

Resumable: a <KO>.fa that already exists is kept (not rewritten), and once a full pass completes a .ko_fastas_complete marker is written so a later call returns the existing files without re-scanning genes_pep (the costly part). force ignores both and rebuilds. The resume assumes the same inputs as the original run (as the per-KO HMM resume in build_hmm_library does); use force after changing organism_gene_ko/organisms.

build_ko_hmm#

Cluster, align and train a profile HMM for one KO's multi-FASTA.

Single-sequence KOs skip CD-HIT/MAFFT (a lone sequence is its own alignment). seq_identity=-1 skips CD-HIT. All (deduplicated) sequences are kept — memory on large KOs is bounded by switching MAFFT to PartTree, not by dropping sequences. fast uses MAFFT FFT-NS-2 (fast progressive) rather than --auto's slow iterative refinement. MAFFT switches to memory-light PartTree once an alignment is predicted to be too memory-heavy: by default from its DP cost (n_seqs × mean_len² — long proteins cost far more than the same residue count in short ones) against a RAM-derived budget (_auto_cost_budget). Passing parttree_residues overrides this with a simple residue-count cutoff. verbose logs (via the logging module, INFO/DEBUG) which tool is running for this KO, sequence counts at each stage, timings, and the tools' own output. Returns out_hmm.

build_reference_model#

build_reference_model(reactions: list[KeggReaction], compounds: list[KeggCompound])

Assemble the gene-free KEGG reference model from parsed records.

Only metabolites actually used by a reaction are added. Reactions carry KEGG annotations (reaction id, KO ids, EC codes, pathways) but no genes/GPRs. Bounds are (-1000, 1000) for reversible reactions and (0, 1000) otherwise.

download_kegg_dump#

download_kegg_dump(dest: str | Path, *, files: tuple[str, ...]=DEFAULT_FILES, base_url: str=BASE_URL, host: str=KEGG_HOST, auth: tuple[str, str] | None=None, netrc_path: str | Path | None=None, force: bool=False, verbose: bool=True, progress: bool=True)

Fetch and arrange a complete KEGG dump into dest.

Convenience wrapper chaining fetch_kegg_files and extract_kegg_dump. Returns the logical-name -> path mapping of the arranged dump, ready for raven_toolbox.reconstruction.kegg.parse_kegg_dump. progress (on by default) shows download/extraction bars; pass progress=False for non-interactive runs.

extract_kegg_dump#

extract_kegg_dump(dest: str | Path, *, progress: bool=False)

Extract and arrange the downloaded archives into the flat dump layout.

Mirrors fetch_keggdb.sh's extract step: untar the *.tar.gz archives, gunzip the *.pep.gz proteomes, lift the needed files out of their sub-directories, and concatenate compound+glycan and the two proteomes. Tar extraction uses the data filter (no path traversal). Returns a mapping of logical name -> path for the files produced. progress shows a bar over the untar step and a byte bar per proteome gunzip.

Network-free, so this is the unit-tested core; download_kegg_dump chains fetch_kegg_files in front of it.

fetch_kegg_files#

fetch_kegg_files(dest: str | Path, *, files: tuple[str, ...]=DEFAULT_FILES, base_url: str=BASE_URL, host: str=KEGG_HOST, auth: tuple[str, str] | None=None, netrc_path: str | Path | None=None, force: bool=False, verbose: bool=True, progress: bool=True)

Download the raw KEGG archives into dest (basenames). Returns the paths.

Existing files are skipped unless force=True (the script's wget -N intent, simplified to skip-if-present). progress shows a per-file byte download bar (sized from the response Content-Length); verbose keeps the fetching/skip INFO log lines.

flag_set#

flag_set(rxn_flags: pd.DataFrame | None, column: str)

Reaction ids whose column flag is truthy (handles bool or TSV strings).

get_kegg_model_for_organism#

get_kegg_model_for_organism(organism_id: str, reference_model: cobra.Model, ko_reaction: pd.DataFrame, organism_gene_ko: pd.DataFrame, *, rxn_flags: pd.DataFrame | None=None, taxonomy: str | Path | None=None, keep_spontaneous: bool=True, keep_undefined_stoich: bool=True, keep_incomplete: bool=True, keep_general: bool=False)

Reconstruct a draft model for a KEGG species from its KO annotations.

Parameters

organism_id: Three/four-letter KEGG organism code (e.g. "eco"), or "eukaryotes"/"prokaryotes" for a whole-domain model (requires taxonomy). Matched case-insensitively.
reference_model: The gene-free KEGG reference model (from build_reference_model).
ko_reaction, organism_gene_ko, rxn_flags: The relational tables from build_kegg_tables (or read back with read_kegg_table).
taxonomy: Path to the KEGG taxonomy file; required only for domain mode.
keep_spontaneous, keep_undefined_stoich, keep_incomplete, keep_general: Quality filters (RAVEN's keep*). A reaction flagged in rxn_flags is dropped unless its keep flag is set; this takes precedence over having genes. Spontaneous reactions are additionally kept without genes when keep_spontaneous is true.

Returns

cobra.Model: A copy of the reference restricted to the organism's reactions, with GPRs built and kegg.genes annotations on the genes.

get_kegg_model_for_organism_from_artefacts#

Load the published 3b.2 artefacts from artefact_dir and build the model.

Reads reference_model.yml.gz and the ko_reaction/organism_gene_ko/ rxn_flags gzipped-TSV tables, then calls get_kegg_model_for_organism. If artefact_dir is None the published artefacts are fetched/cached via raven_toolbox.data.ensure_kegg_data (version selects the release).

For a whole-domain model (organism_id = "prokaryotes"/"eukaryotes") the KEGG taxonomy artefact is resolved automatically — from artefact_dir if present, otherwise fetched via raven_toolbox.data.ensure_kegg_taxonomy (it is a separate artefact, not part of the core set) — unless an explicit taxonomy path is given.

get_kegg_model_from_sequences#

get_kegg_model_from_sequences(fasta: str | Path, reference_model: cobra.Model, ko_reaction: pd.DataFrame, library: str | Path, *, rxn_flags: pd.DataFrame | None=None, model_id: str | None=None, cutoff: float=1e-30, min_score_ratio_ko: float=0.3, min_score_ratio_g: float=0.9, keep_spontaneous: bool=True, keep_undefined_stoich: bool=True, keep_incomplete: bool=True, keep_general: bool=False, threads: int=max(1, (os.cpu_count() or 2) - 1), hmmsearch: str | Path | None=None)

Reconstruct a draft model for a proteome by HMM-searching the KO library.

Searches the library against fasta (3b.3), assigns KOs (assign_kos), and assembles the model against reference_model / ko_reaction. Genes are the query proteome's identifiers.

get_kegg_model_from_sequences_with_artefacts#

Load reference model + tables from artefact_dir and run the HMM query.

If artefact_dir / library are None they are fetched/cached via raven_toolbox.data.ensure_kegg_data / raven_toolbox.data.ensure_kegg_hmm_library (domain selects the prok/euk library; version the release).

get_model_from_homology#

get_model_from_homology(models, hits: pd.DataFrame, model_for: str, *, preferred_order=None, strictness: int=1, map_direction: str='new_to_old', complex_policy: str='flag', only_genes_in_models: bool=False, max_evalue: float=1e-30, min_align_len: int=100, min_identity: float=40, review_identity: float | None=None, keep_gene_free: bool=False, preserve_notes: bool=False)

Build a draft model for model_for by transferring reactions from templates.

Other parameters that materially change the result:

strictness (default 1) picks the hit-filtering policy:
- 1: bidirectional — a template gene and a new-organism gene must have an accepted hit in both directions, rather than in just one.
- 2: one-directional — trust a hit found in a single direction, picked by map_direction.
- 3: bidirectional and best-hits-only — before matching, each gene keeps only its best hits (by bitscore) in each direction, all of them when several tie, then the bidirectional rule is applied.
map_direction picks which one-directional search to trust when strictness is 2: "new_to_old" (default, hits found searching from model_for) or "old_to_new".
complex_policy decides what happens to an AND-linked subunit with no ortholog: "flag" (default) keeps the reaction with a placeholder gene for later curator review, "keep" drops just that subunit, "drop" drops the whole reaction.
only_genes_in_models restricts the hits table to hits whose template gene appears in the given models before mapping orthologs, so a gene the templates lack cannot be a best hit. It raises ValueError when no hit of a search direction between model_for and a template passes the cutoffs and this filter.
keep_gene_free (default False) also transfers the template reactions that have no GPR (spontaneous, transport, exchange, artificial). By default they are dropped, as in RAVEN, because nothing supports them in the new organism. Set it when the template is a curated reference whose gene-free reactions belong in every derived model.
preserve_notes (default False) keeps each transferred reaction's own notes (references, confidence score, free text) and only adds homology_source. By default they are replaced by a note that the reaction was included by this function.
preferred_order breaks ties when more than one template maps the same new-organism gene: the earlier template in the list wins that gene (default: the order of models).

Organism ids in hits are case-sensitive. An id that matches model_for or a template id only when case is ignored raises ValueError; hits involving other organisms are ignored.

Defaults for the three filters come from a calibration against KEGG orthology <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/homology-cutoff-calibration.md>_ on raven-docs, scored with precision weighted above recall (a wrongly transferred reaction is harder to undo than a missing one):

min_identity 40 is the binding filter and the measured optimum.
min_align_len 100 replaces RAVEN's 200, which discarded real orthologs for no gain in precision. Anything at or below 150 measured the same; 200 cost 1-2 points on every organism tested. A hit passes when its e-value is below max_evalue and its alignment length and identity are at least min_align_len and min_identity, as in RAVEN.
max_evalue makes no difference anywhere between 1e-4 and 1e-50 -- the other two filters have already excluded whatever it would exclude -- so 1e-30 is kept for continuity with RAVEN rather than for any measured effect.

review_identity (e.g. 25) collects what those thresholds turn away. Any reaction that would have transferred at the looser identity is reported in HomologyResult.candidates with its strongest supporting hit, and is not added to the model. The two error types are not symmetric -- a missing reaction can be gap-filled, a wrong one is hard to find and harder to remove -- so the filters stay strict, but the evidence for the near misses is handed to the curator instead of being discarded silently.

hits_from_supporth#

Build a homology hits table from a SuppOrth suppOrth.tsv.

Only the pairs between organism_id and another species are kept. Each pair is written in both directions, so every strictness of raven_toolbox.reconstruction.homology.get_model_from_homology accepts it, and strictness 1 and 2 give the same result.

SuppOrth decides orthology by agreement between predictors, so the filters that apply to it are in this function, and the hit metrics are placeholders that pass the alignment filters of get_model_from_homology (evalue 0, identity 100, align_len 1000, ppos 100). bitscore is SuppOrth's own value (0 when it has none), which strictness=3 ranks by. Because the identity is a placeholder, review_identity reports no candidates for these hits.

Parameters

source: Path to suppOrth.tsv, or the table as a DataFrame. A DataFrame read with pandas defaults has protein ids such as NA turned into NaN; read it with keep_default_na=False or pass the path.
organism_id: SuppOrth species label of the organism to build a model for: the FASTA file name without extension, or the --labels value. It becomes from_id / to_id unless species_to_model maps it to another id.
species_to_model: SuppOrth species label to model id. A label mapped to None is left out. When given, every label in the table must be a key. When None, labels are used as model ids.
min_support: Minimum number of predictors that found the pair (Lv_support). SuppOrth has already filtered the table with its own --supportLv and --identity_filter; this filter is applied in addition.
require_tools: Predictors that must all have found the pair, by column name ("orthofinder", "sonicparanoid", "broccoli", "oma").
min_identity, max_evalue: Filters on SuppOrth's alignment values. A pair without an alignment hit fails both.
protein_to_gene: Protein id to model gene id, as a dict, a JSON file, or a two-column tab-separated file (lines starting with # are skipped). Proteins missing from the map keep their id, and their number is reported in a warning. A gene pair takes the highest support over its protein pairs, the union of their predictors, the highest identity and bitscore, and the lowest e-value; these values can come from different protein pairs.

Returns

pandas.DataFrame: HIT_COLUMNS followed by support, supported_by, identity_raw, evalue_raw and bitscore_raw (SuppOrth's values). attrs["source"] is "SuppOrth".

HomologyResult#

Result of get_model_from_homology.

Attributes

model: The draft cobra.Model.
gene_map: {model_id: {template_gene: [new_gene, ...]}} ortholog mapping used.
candidates: Reactions that just missed the identity threshold, with the best hit supporting each, sorted strongest first. None unless review_identity was given. These are not in model: they are the near misses, for a curator to accept or reject deliberately.

KeggCompound#

A metabolite parsed from the KEGG compound flat file.

KeggKO#

A KEGG Orthology entry: its name and the organism genes assigned to it.

KeggReaction#

A reaction parsed from the KEGG reaction flat file.

make_ortholog_hits#

make_ortholog_hits(ortholog_pairs: Iterable[tuple[str, str]], source_model_id: str, target_id: str)

Build a bidirectional hits table from a predefined ortholog list.

Each (source_gene, target_gene) pair is emitted in both directions with sentinel metrics (evalue 0, identity 100, align_len 1000, bitscore 1000, ppos 100) so every pair passes any reasonable filter. Lets a known ortholog mapping feed get_model_from_homology with no BLAST run — also the testing entry point.

Parameters

ortholog_pairs: Iterable of (source_gene, target_gene) — source = template/model organism, target = the organism being built.
source_model_id: ID of the template model the source genes belong to.
target_id: ID of the organism to build a model for (model_for).

organism_domains#

organism_domains(path: str | Path)

Return {organism_id: domain} (the outermost category).

organisms_in_domain#

organisms_in_domain(path: str | Path, domain: str)

Organism ids whose outermost category matches domain (case-insensitive).

Accepts a prefix, so "prok" matches "Prokaryotes".

parse_hmmsearch_tblout#

parse_hmmsearch_tblout(text: str)

Parse hmmsearch --tblout text into a [ko, gene, evalue] table.

With the profile library as the query, column 1 (target name) is the proteome gene, column 3 (query name) is the KO, and column 5 is the full-sequence E-value.

parse_kegg_compounds#

parse_kegg_compounds(kegg_dir: str | Path)

Parse <kegg_dir>/compound (+ optional compound.inchi) into records.

parse_kegg_dump#

parse_kegg_dump(kegg_dir: str | Path, out_dir: str | Path, *, version: str | None=None, progress: bool=False)

Parse a full KEGG dump into the reference model + tables and write them out.

Writes reference_model.yml.gz (gzipped RAVEN/cobra YAML) plus the gzipped-TSV tables into out_dir and returns {name: path} for everything written. The large organism_gene_ko table is streamed to disk (see stream_organism_gene_ko) rather than built in memory, so this scales to the full KEGG database; the small derived tables are built in memory.

When version is given (e.g. "kegg116") the output filenames are version-prefixed (<version>_<name>), matching the published release assets; the returned dict keys stay the logical table names.

progress reports each stage and shows a byte-progress bar over the dominant organism_gene_ko (ko) streaming pass.

parse_kegg_kos#

parse_kegg_kos(kegg_dir: str | Path, *, keep: set[str] | None=None)

Parse <kegg_dir>/ko into KeggKO records (name + organism genes).

keep limits parsing to those KO ids (e.g. only KOs linked to reactions), mirroring RAVEN's koList argument — the gene lists are huge, so this is the usual call.

parse_kegg_reactions#

parse_kegg_reactions(kegg_dir: str | Path)

Parse <kegg_dir>/reaction into KeggReaction records.

Reversibility is taken from the equation arrow and, when reaction_mapformula.lst is present, refined to mark reactions that are irreversible across all KEGG maps (see _irreversible_from_mapformula).

parse_taxonomy#

parse_taxonomy(path: str | Path)

Return {organism_id: [category, ...]} from outermost to innermost.

parse_taxonomy_records#

parse_taxonomy_records(path: str | Path)

Return [(organism_id, name, [category, ...]), ...] in file order.

name is the organism's scientific name with any trailing KEGG parenthetical kept (e.g. "Homo sapiens (human)"), matching RAVEN getPhylDist; the category list is the lineage from outermost (domain) to innermost. Order is preserved so it aligns with the rows/columns of phyl_dist.

phyl_dist#

phyl_dist(path: str | Path, *, only_in_kingdom: bool=False)

Compute the KEGG phylogenetic distance matrix from the taxonomy file.

Faithful port of RAVEN getPhylDist. For organisms i and j with category lineages of lengths Li, Lj, the distance is (Li - Lj) + min(Li, Lj) - k where k is the deepest level (within the shorter lineage, comparing position by position from the root) at which the two categories coincide — 1 if none coincide. This reproduces RAVEN's exact values (including their asymmetry and occasional negatives) so the output matches MATLAB keggPhylDist that GECKO was built against.

only_in_kingdom sets the distance to +inf between organisms in different domains (RAVEN's onlyInKingdom variant); the matrix is then returned as float.

Cost is O(N²) in time and memory (N = number of KEGG organisms, ~10⁴), so this is a maintainer/GECKO-side generation step; persist the result rather than rebuilding.

PhylDist#

KEGG phylogenetic distance structure — RAVEN getPhylDist / keggPhylDist.

ids/names are the KEGG organism ids and scientific names (names keep RAVEN's trailing parenthetical; consumers that need them cleaned, e.g. geckopy, strip it). dist_matrix[i, j] is RAVEN's distance from ids[i] to ids[j]. It is a faithful reproduction of RAVEN's metric, which is asymmetric and may be negative; GECKO consumes it only by taking the closest (minimum-distance) organism, so the exact values matter for parity with MATLAB but the sign/asymmetry are harmless there.

read_kegg_table#

read_kegg_table(path: str | Path)

Read a KEGG table written by write_kegg_tables or

stream_organism_gene_ko.

Compression is inferred from the suffix; all published tables are gzipped TSV (.tsv.gz), and a version-prefixed name (kegg116_<name>.tsv.gz) reads just the same.

run_blast#

Bidirectional BLAST+ between an organism and template organisms.

model_ids and ref_fastas are parallel: ref_fastas[i] is the reference proteome for template model_ids[i].

Returns the hits DataFrame (filtered at evalue). Requires BLAST+ (blastp, makeblastdb).

The default matches MATLAB RAVEN's getBlast, which hardcodes -evalue 10e-5 -- that is 1e-4, not 1e-5. Unlike RAVEN's, this one is an argument, so a stricter search stays available; it's also a loose pre-filter only, dominated downstream by raven_toolbox.reconstruction.homology.homology.get_model_from_homology's much stricter max_evalue/min_align_len/min_identity.

run_diamond#

run_diamond(organism_id: str, fasta: str | Path, model_ids: Sequence[str], ref_fastas: Sequence[str | Path], *, evalue: float=0.001, threads: int=max(1, (os.cpu_count() or 2) - 1), sensitivity: str='--more-sensitive', diamond: str | Path | None=None)

Bidirectional DIAMOND between an organism and template organisms.

model_ids and ref_fastas are parallel: ref_fastas[i] is the reference proteome for template model_ids[i]. sensitivity is passed straight through as a DIAMOND CLI flag (e.g. "--more-sensitive", "--sensitive", or "" to use DIAMOND's default fast mode).

Returns the hits DataFrame. Requires DIAMOND.

The default matches MATLAB RAVEN's getDiamond, which passes no --evalue at all and so inherits DIAMOND's own default of 1e-3. Stated explicitly here rather than left implicit, since the two aligners disagree: run_blast uses 1e-4, following getBlast. Neither cutoff has been calibrated against data on either side -- unlike the KO-assignment cutoffs, where measurement did justify diverging from RAVEN.

run_hmmsearch#

run_hmmsearch(fasta: str | Path, library: str | Path, *, threads: int=max(1, (os.cpu_count() or 2) - 1), hmmsearch: str | Path | None=None)

Search the profile library against proteome fasta; return tblout text.

One hmmsearch of the whole concatenated multi-profile library (query) against fasta (target) — the fast search direction, parallelised with --cpu, and no hmmpress needed. -Z is fixed to the profile count so the per-hit E-values match the convention assign_kos is calibrated against (identical to a hmmscan against the same library).

run_supporth#

run_supporth(organism_id: str, fasta: str | Path, model_ids: Sequence[str], ref_fastas: Sequence[str | Path], outdir: str | Path, *, tools: str | Sequence[str]='all', support_lv: int | None=None, identity_filter: float=80.0, alignment: str='diamond', threads: int=max(1, (os.cpu_count() or 2) - 1), extra_args: Sequence[str]=(), supporth: str | os.PathLike | None=None, **reader_kwargs)

Run SuppOrth on an organism and template proteomes and read its consensus.

model_ids and ref_fastas are parallel: ref_fastas[i] is the proteome of template model_ids[i], with FASTA ids equal to that model's gene ids (or mapped to them with protein_to_gene). The ids are passed to SuppOrth as --labels, so the table uses them as species labels.

support_lv and identity_filter are SuppOrth's --supportLv and --identity_filter: a pair is kept if at least support_lv predictors found it, or if its identity is above identity_filter. support_lv defaults to 3, or to the number of tools when fewer are selected. extra_args are appended to the command, for example ["--resume"] or ["--oma", "path/to/Output"]. Output is written to outdir.

SuppOrth requires protein ids to be unique over all proteomes. When two FASTA files share an id (for example two templates of the same species), every id is given a <label>__ prefix in copies written to outdir/prefixed_inputs, and the prefix is removed from the returned hits. outdir/suppOrth.tsv then contains the prefixed ids. SuppOrth's console output is written to outdir/run_supporth.log.

Remaining keyword arguments go to hits_from_supporth.

Requires the suppOrth command, found from supporth, then the RAVEN_SUPPORTH environment variable, then PATH. SuppOrth's predictors do not run on Windows; use WSL or a Linux host, and read the result with hits_from_supporth.

stream_organism_gene_ko#

stream_organism_gene_ko(kegg_dir: str | Path, keep: set[str], ogk_path: str | Path, *, chunk_rows: int=1000000, progress: bool=False)

Stream the ko file to a sorted, gzipped organism_gene_ko.tsv.gz.

keep restricts output to that set of KO ids (e.g. only KOs linked to reactions), mirroring keep on parse_kegg_kos.

Real KEGG has ~9M gene↔KO associations — far too many to hold in memory as a DataFrame. Rows are sorted by (organism, gene) before writing: gene IDs from one organism share long common prefixes (locus tags, numeric runs), so sorting makes them adjacent and helps the compressor; the order also matches the by-organism query pattern in get_kegg_model_for_organism. Gzip (not xz) keeps the table readable by MATLAB's built-in gunzip with no external tool, at a modestly larger size.

The sort is an external merge sort bounded to chunk_rows rows in memory at a time (sorted runs spooled to gzipped temp files, then merged with heapq.merge), so peak memory stays flat regardless of KEGG size. Only the small ko_names table (one row per KO) is held in full and returned.

progress shows a tqdm byte-progress bar over the (large) ko read pass.

validate_hits#

validate_hits(hits: pd.DataFrame)

Check a hits DataFrame has the required columns; return it unchanged.

write_kegg_tables#

write_kegg_tables(tables: dict[str, pd.DataFrame], out_dir: str | Path, *, prefix: str='')

Write each table as a gzipped TSV (<prefix><name>.tsv.gz) into out_dir.

Gzipped TSV is the dependency-free cross-language format shared with MATLAB RAVEN (see docs/kegg_data_format.md) — readable by MATLAB's built-in gunzip with no external tool. prefix version-tags the filenames (e.g. kegg116_). Returns the written paths.