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 from a CSV path or a DataFrame (not Excel).

build_hmm_library

Build the HMM library of a domain ("prokaryotes" or "eukaryotes").

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, named by their basenames, and return the paths.

flag_set

Reaction ids whose column flag is truthy.

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 reference 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 the reference model and 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 and the 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 and 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 of RAVEN getPhylDist and keggPhylDist.

read_kegg_table

Read a KEGG table from write_kegg_tables or stream_organism_gene_ko.

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 that a hits DataFrame has the required columns and 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 and 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 contributed a gene, as in the HMM branch of RAVEN getKEGGModelForOrganism. The organism-annotation path keeps the full reference KO list, and the default is False.

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 by gene E-value matrix. An E-value of 0 is replaced by 1e-250. The steps are:

  1. keep hits with evalue <= cutoff;

  2. min_score_ratio_ko: within a KO, drop genes whose log(evalue)/log(best_evalue_in_KO) < min_score_ratio_ko (removes weak members of a KO);

  3. min_score_ratio_g: within a gene, drop KOs whose log(evalue)/log(best_evalue_for_gene) < min_score_ratio_g (removes the weaker KOs of a gene that has a better KO).

A smaller E-value is a better hit. All kept E-values are below 1, so their logs are negative, the best (smallest) hit has ratio 1, and weaker hits have a smaller positive ratio.

The defaults were calibrated by cross-validation against the KEGG gene and KO annotation of four organisms that span the prokaryote and eukaryote libraries and well- and less-studied species (S. cerevisiae, Cyanidioschyzon merolae, E. coli, Mycoplasma genitalium). Annotated genes have a median E-value of 1e-100 to 1e-155, and spurious hits cluster at about 1e-8, a separation of ~20 orders of magnitude. RAVEN's cut-off of 1e-50 lies within the E-value range of the annotated genes and drops real but divergent hits. For the divergent minimal genome (M. genitalium, mge) the gene-to-KO recall is 0.95 at 1e-30 and 0.84 at 1e-50, with no gain in rejecting noise (noise E-values are far weaker). The default cut-off is therefore 1e-30, which recovers those hits and is ~22 orders of magnitude above the noise level. min_score_ratio_g (default 0.9) sets the precision, because it resolves genes with several KOs. Against RAVEN's settings (1e-50 and 0.8), the recall for mge rises from 0.84 to 0.94. min_score_ratio_ko values of 0.0, 0.3 and 0.5 changed precision and recall by at most 0.02 across the four organisms; 0.3 is kept as in RAVEN.

Parameters

hits : pandas.DataFrame

Hit table with the columns ko, gene and evalue, as returned by parse_hmmsearch_tblout. For repeated (ko, gene) pairs the smallest E-value is used.

cutoff : float

Largest E-value that is kept. Defaults to 1e-30. It must be below 1; a value of 1 or more raises a ValueError.

min_score_ratio_ko : float

Smallest ratio log(evalue)/log(best_evalue_in_KO) that a gene must reach to stay in a KO. Defaults to 0.3. A value of 0 keeps all genes of a KO.

min_score_ratio_g : float

Smallest ratio log(evalue)/log(best_evalue_for_gene) that a KO must reach to stay assigned to a gene. Defaults to 0.9.

blast_from_table#

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

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

Must contain the HIT_COLUMNS columns.

Parameters

source : str or Path or pandas.DataFrame

A DataFrame, or the path of a CSV file. The columns from_gene and to_gene of a CSV are read as text, and numbers are parsed exactly. A table that lacks a column of HIT_COLUMNS raises ValueError.

build_hmm_library#

build_hmm_library(organism_gene_ko: pd.DataFrame, genes_pep: str | Path, taxonomy: str | Path, out_dir: str | Path, *, domain: str, seq_identity: float=0.9, parttree_residues: int | None=None, threads: int | None=None, fast: bool=True, verbose: bool=False, progress: bool=False, force: bool=False, concatenate: bool=True, cdhit: str | Path | None=None, mafft: str | Path | None=None, hmmbuild: str | Path | None=None)

Build the HMM library of a domain ("prokaryotes" or "eukaryotes").

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, which hmmsearch searches directly without hmmpress. Returns {"hmms": [...], "library": path | None}.

The function is run by the maintainer once per KEGG release and needs CD-HIT, MAFFT and hmmbuild. It skips KOs that already have an .hmm and reuses an existing multi-FASTA (see build_ko_fastas), so a rerun after a failure continues from the existing files. force rebuilds everything.

Parameters

organism_gene_ko : pandas.DataFrame

Gene-to-KO table with the columns organism, gene and ko. See build_ko_fastas.

genes_pep : str or Path

Multi-FASTA file of the KEGG proteomes (genes.pep).

taxonomy : str or Path

KEGG taxonomy file, from which the organism codes of domain are taken.

out_dir : str or Path

Output directory. It receives the sub-directories fasta and hmms and, if concatenate, library.hmm.

domain : str

Domain whose organisms are used, "prokaryotes" or "eukaryotes". A ValueError is raised if the taxonomy has no organism for it.

seq_identity : float

CD-HIT identity threshold, in (0.4, 1.0]. Defaults to 0.9. A value of -1 skips CD-HIT. Passed to build_ko_hmm.

parttree_residues : int or None

Cutoff on the residue count above which MAFFT uses PartTree. When None (default), the DP-cost criterion of build_ko_hmm is used.

threads : int or None

Number of threads for CD-HIT, MAFFT and hmmbuild, passed to build_ko_hmm. When None (default), every CPU the process may use.

fast : bool

If True (default), MAFFT uses FFT-NS-2. If False, it uses --auto. Passed to build_ko_hmm.

verbose : bool

If True, log the per-KO steps at INFO and DEBUG. Defaults to False.

progress : bool

If True, show tqdm bars for the multi-FASTA and per-KO HMM passes, with an N of M counter for the per-KO pass. It is independent of verbose. With both on, per-KO log lines are routed through tqdm.write so that they do not overwrite the bar. Defaults to False.

force : bool

If True, rebuild every multi-FASTA and every .hmm. If False (default), existing files are reused.

concatenate : bool

If True (default), write the per-KO HMMs, in sorted order, to out_dir/library.hmm. If False, or if no HMM was built, library is None in the result.

cdhit : str or Path or None

Path of the cd-hit executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary.

mafft : str or Path or None

Path of the mafft executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary.

hmmbuild : str or Path or None

Path of the hmmbuild executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary.

build_kegg_tables#

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

Build the minimal relational tables from parsed records.

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

Parameters

reactions : list of KeggReaction

Records from parse_kegg_reactions. They give the tables ko_reaction and rxn_flags.

kos : list of KeggKO

Records from parse_kegg_kos. They give the tables ko_names and organism_gene_ko.

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.

The sequences are located through an index of byte offsets built in one pass over genes_pep. KOs without a sequence in genes_pep are skipped and get no file. Returns {ko: path} for the files written.

A <KO>.fa that already exists is kept and not rewritten. Once a full pass completes, a .ko_fastas_complete marker is written, and a later call returns the existing files without scanning genes_pep again. force ignores both and rebuilds. The resume assumes the same inputs as the original run, as does the per-KO HMM resume in build_hmm_library; force is needed after organism_gene_ko or organisms change.

Parameters

organism_gene_ko : pandas.DataFrame

Gene-to-KO table with the columns organism, gene and ko. The sequence id of a row is organism:gene, as in the genes_pep headers.

genes_pep : str or Path

Multi-FASTA file of the KEGG proteomes (genes.pep). The id of a record is the first whitespace-delimited token of its > header.

out_dir : str or Path

Directory for the <KO>.fa files. It is created, with parents, if missing.

organisms : set of str or None

Organism codes to keep, for example the codes of one domain for the prokaryote and eukaryote libraries. When None (default), all rows of organism_gene_ko are used.

progress : bool

If True, show a tqdm bar over the per-KO writing pass. Defaults to False.

force : bool

If True, rebuild every <KO>.fa and ignore the completion marker. If False (default), existing files and the marker are reused.

build_ko_hmm#

build_ko_hmm(ko_fasta: str | Path, out_hmm: str | Path, *, seq_identity: float=0.9, parttree_residues: int | None=None, threads: int | None=None, fast: bool=True, verbose: bool=False, cdhit: str | Path | None=None, mafft: str | Path | None=None, hmmbuild: str | Path | None=None)

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

A KO with one sequence skips CD-HIT and MAFFT, because a single sequence is its own alignment; MAFFT is also skipped when one sequence remains after CD-HIT. CD-HIT removes near-identical sequences at seq_identity, and no other sequences are dropped. MAFFT switches to PartTree, which needs less memory, once an alignment is predicted to need too much: by default from its DP cost (n_seqs × mean_len²) against a RAM-derived budget (_auto_cost_budget). The cost of an alignment of long proteins is higher than that of the same residue count in short proteins. Returns out_hmm.

Parameters

ko_fasta : str or Path

Multi-FASTA file of the sequences of one KO. A ValueError is raised if it contains no sequence.

out_hmm : str or Path

Path of the profile HMM to write. The parent directory is created if missing, and the file stem names the profile.

seq_identity : float

CD-HIT identity threshold, in (0.4, 1.0]. Defaults to 0.9. A value of -1 skips CD-HIT. A value outside the range raises a ValueError when CD-HIT runs, that is for a KO with more than one sequence.

parttree_residues : int or None

Cutoff on the residue count of the clustered sequences. An alignment with more residues uses MAFFT PartTree. When None (default), the DP-cost criterion above is used instead.

threads : int or None

Number of threads for CD-HIT, MAFFT and hmmbuild. When None (default), every CPU the process may use.

fast : bool

If True (default), MAFFT uses FFT-NS-2 (--retree 2 --maxiterate 0), a fast progressive alignment. If False, it uses --auto. PartTree takes precedence over both.

verbose : bool

If True, log via the logging module at INFO and DEBUG which tool is running for this KO, the sequence counts at each stage, the timings and the output of the tools. Defaults to False.

cdhit : str or Path or None

Path of the cd-hit executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary. It is resolved only when CD-HIT runs.

mafft : str or Path or None

Path of the mafft executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary. It is resolved only when MAFFT runs.

hmmbuild : str or Path or None

Path of the hmmbuild executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary.

build_reference_model#

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

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

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

Parameters

reactions : list of KeggReaction

Records from parse_kegg_reactions. A record whose cached stoichiometry is empty is parsed again from its equation.

compounds : list of KeggCompound

Records from parse_kegg_compounds. They give the name, formula and annotations of the metabolites. A metabolite without a record is named by its id.

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.

Runs fetch_kegg_files and then extract_kegg_dump. Returns the mapping of logical name to path of the arranged dump, the input of raven_toolbox.reconstruction.kegg.parse_kegg_dump.

Parameters

dest : str or Path

Directory that receives the dump. It is created, with parents, if missing.

files : tuple of str

KEGG FTP paths to download, relative to base_url. Defaults to DEFAULT_FILES. See fetch_kegg_files.

base_url : str

Base URL of the server. Defaults to "https://ftp.kegg.net".

host : str

Host name looked up in the .netrc file. Defaults to "ftp.kegg.net". Unused when auth is given.

auth : tuple of (str, str) or None

(user, password) for the server. When None (default), the credentials are read from the .netrc file.

netrc_path : str or Path or None

Path of the .netrc file. When None (default), ~/.netrc is used. Ignored when auth is given.

force : bool

If True, archives that already exist in dest are downloaded again. If False (default), they are skipped.

verbose : bool

If True (default), log the download, skip and extraction steps at INFO level.

progress : bool

If True (default), show download and extraction progress bars.

Returns

dict of str to Path

Logical name to path of each arranged file that exists in dest.

extract_kegg_dump#

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

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

Mirrors the extract step of fetch_keggdb.sh: untars the *.tar.gz archives, gunzips the *.pep.gz proteomes, moves the needed files out of their sub-directories, and concatenates compound and glycan, and the two proteomes. Tar extraction uses the data filter, which rejects paths outside dest. Returns a mapping of logical name to path for the files produced. progress shows a bar over the untar step and a byte bar per proteome gunzip.

The function makes no network access; download_kegg_dump runs fetch_kegg_files before it.

Parameters

dest : str or Path

Directory that holds the downloaded archives. The extracted files are written into it, and the archives and sub-directories are deleted. A FileNotFoundError is raised if the reaction, compound or ko archive did not yield its file.

progress : bool

If True, show a progress bar over the untar step and a byte bar for each proteome gunzip. Defaults to False.

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, named by their basenames, and return the paths.

Existing files are skipped unless force=True. progress shows a per-file byte download bar, sized from the response Content-Length; verbose enables the fetching and skip INFO log lines.

Parameters

dest : str or Path

Directory that receives the archives. It is created, with parents, if missing.

files : tuple of str

KEGG FTP paths to download, relative to base_url. Each file is saved under its basename. Defaults to DEFAULT_FILES: the reaction, compound, glycan and ko archives, the eukaryote and prokaryote proteomes, and the taxonomy file.

base_url : str

Base URL of the server. Defaults to "https://ftp.kegg.net". A trailing slash is removed.

host : str

Host name looked up in the .netrc file. Defaults to "ftp.kegg.net". It is unused when auth is given.

auth : tuple of (str, str) or None

(user, password) for HTTP basic or digest authentication. When None (default), the credentials of host are read from the .netrc file. A FileNotFoundError is raised if that file does not exist, and a ValueError if it cannot be parsed or has no complete entry for host.

netrc_path : str or Path or None

Path of the .netrc file. When None (default), ~/.netrc is used. Ignored when auth is given.

force : bool

If True, files that already exist in dest are downloaded again. If False (default), they are skipped.

verbose : bool

If True (default), log a fetching line for each download and a skip line for each existing file at INFO level.

progress : bool

If True (default), show a byte progress bar for each download. The total comes from the Content-Length header; without it the bar shows a byte count.

Returns

list of Path

Paths of the files in dest, in the order of files, including skipped files.

flag_set#

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

Reaction ids whose column flag is truthy.

A flag is truthy when it is True or when its string form is "true" or "1" in any case.

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- or four-letter KEGG organism code (for example "eco"), or "eukaryotes" or "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 in 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 without genes are also kept 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#

get_kegg_model_for_organism_from_artefacts(organism_id: str, artefact_dir: str | Path | None=None, *, version: str | None=None, taxonomy: str | Path | None=None, **kwargs)

Load the published reference artefacts from artefact_dir and build the model.

Reads reference_model.yml.gz and the gzipped TSV tables ko_reaction, organism_gene_ko and rxn_flags, then calls get_kegg_model_for_organism. If artefact_dir is None, the published artefacts are fetched and cached via raven_toolbox.data.ensure_kegg_data (version selects the release). For a whole-domain model (organism_id "prokaryotes" or "eukaryotes") the KEGG taxonomy artefact is taken from artefact_dir if present, otherwise fetched via raven_toolbox.data.ensure_kegg_taxonomy, because it is a separate artefact outside the core set. An explicit taxonomy path takes precedence.

Parameters

organism_id : str

Three- or four-letter KEGG organism code (for example "eco"), or "eukaryotes" or "prokaryotes" for a whole-domain model. Matched case-insensitively.

artefact_dir : str or Path or None

Directory with the artefact files reference_model.yml.gz, ko_reaction.tsv.gz, organism_gene_ko.tsv.gz and rxn_flags.tsv.gz, each with or without a <version>_ prefix. When None (default), the published artefacts are fetched and cached.

version : str or None

KEGG release to fetch. Used only when artefact_dir is None, or when the taxonomy of a domain model is fetched. When None (default), the default release of raven_toolbox.data.ensure_kegg_data is used.

taxonomy : str or Path or None

Path of the KEGG taxonomy file for a domain model. When None (default), it is taken from artefact_dir or fetched. Unused for a single organism.

**kwargs : dict

Keyword arguments forwarded to get_kegg_model_for_organism: the quality filters keep_spontaneous, keep_undefined_stoich, keep_incomplete and keep_general. The arguments rxn_flags and taxonomy are set by this function.

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 | None=None, hmmsearch: str | Path | None=None)

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

Searches the library against fasta, assigns KOs (assign_kos), and assembles the model against reference_model and ko_reaction. Genes are the identifiers of the query proteome. model_id defaults to the stem of the fasta file name.

Parameters

fasta : str or Path

Proteome FASTA file of the query organism.

reference_model : cobra.Model

Gene-free KEGG reference model (from build_reference_model).

ko_reaction : pandas.DataFrame

Table with the columns ko and reaction that maps KOs to reactions.

library : str or Path

Concatenated KO profile-HMM library, searched by run_hmmsearch.

rxn_flags : pandas.DataFrame or None

Per-reaction quality-flag table with a reaction column and any of the boolean columns spontaneous, undefined_stoich, incomplete and general. When None (default), no reaction is flagged.

model_id : str or None

Id of the model. When None (default), the stem of the fasta file name.

cutoff : float

Largest E-value that is kept, passed to assign_kos. Defaults to 1e-30.

min_score_ratio_ko : float

Minimum score ratio of a gene within a KO, passed to assign_kos. Defaults to 0.3.

min_score_ratio_g : float

Minimum score ratio of a KO within a gene, passed to assign_kos. Defaults to 0.9.

keep_spontaneous : bool

If True (default), reactions flagged spontaneous in rxn_flags are kept without genes.

keep_undefined_stoich : bool

If True (default), reactions flagged undefined_stoich in rxn_flags are kept. If False, they are dropped even if they have genes.

keep_incomplete : bool

If True (default), reactions flagged incomplete in rxn_flags are kept. If False, they are dropped even if they have genes.

keep_general : bool

If True, reactions flagged general in rxn_flags are kept. If False (default), they are dropped even if they have genes.

threads : int or None

Number of hmmsearch threads, passed to run_hmmsearch. When None (default), every CPU the process may use.

hmmsearch : str or Path or None

Path of the hmmsearch executable, passed to run_hmmsearch. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary.

get_kegg_model_from_sequences_with_artefacts#

get_kegg_model_from_sequences_with_artefacts(fasta: str | Path, artefact_dir: str | Path | None=None, library: str | Path | None=None, *, domain: str='prokaryotes', version: str | None=None, **kwargs)

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

If artefact_dir or library is None, it is fetched and cached via raven_toolbox.data.ensure_kegg_data or raven_toolbox.data.ensure_kegg_hmm_library. domain selects the prokaryote or eukaryote library and version the release.

Parameters

fasta : str or Path

Proteome FASTA file of the query organism.

artefact_dir : str or Path or None

Directory with the artefact files reference_model.yml.gz, ko_reaction.tsv.gz and rxn_flags.tsv.gz, each with or without a <version>_ prefix. When None (default), the published artefacts are fetched and cached.

library : str or Path or None

Concatenated KO profile-HMM library. When None (default), the library of domain is fetched and cached.

domain : str

Library to fetch when library is None, "prokaryotes" (default) or "eukaryotes". Unused when library is given.

version : str or None

KEGG release to fetch. Used only for the files that are fetched. When None (default), the default release of the data registry is used.

**kwargs : dict

Keyword arguments forwarded to get_kegg_model_from_sequences: model_id, cutoff, min_score_ratio_ko, min_score_ratio_g, keep_spontaneous, keep_undefined_stoich, keep_incomplete, keep_general, threads and hmmsearch. The argument rxn_flags is set by this function.

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.

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.

Defaults for the three filters come from a calibration against KEGG annotations and OMA assignments <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).

Parameters

models : cobra.Model or Sequence[cobra.Model]

Template models. Each template is identified by its model.id, which must match the from_id or to_id values in hits. A single model is accepted as a list of one. Templates are merged in the order they are given, whatever the preferred order is.

hits : pandas.DataFrame

Bidirectional hits table with the columns in HIT_COLUMNS; a missing column raises ValueError. 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. RAVEN raises an error for hits from an organism that no template has as model.id, and matches the template ids case-insensitively.

model_for : str

Id of the organism the draft model is built for, as it appears in hits.

preferred_order : Sequence[str] or None

Breaks ties when more than one template maps the same new-organism gene: the gene is mapped to the earlier template in the list, and a gene of a later template that maps to any gene an earlier template already maps is removed from that template altogether, with its other mappings. A reaction of the later template that cannot work without the removed gene is not transferred. Templates are merged in the order they were given in models, whatever the preferred order is, and a template that preferred_order leaves out contributes nothing. The tie-break applies only when more than one model is given. Default None, which uses the order of models. RAVEN's preferredOrder decides which template receives a gene in the same way and merges the models in the order of models.

strictness : int

Hit-filtering policy. Default 1. Any other value than 1, 2 or 3 raises ValueError.

  • 1: bidirectional. A template gene and a new-organism gene must have an accepted hit in both directions.

  • 2: one-directional. A hit found in a single direction is accepted; the direction is selected 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 : str

Direction of the one-directional search when strictness is 2: "new_to_old" (hits found searching from model_for) or "old_to_new". Default "new_to_old". Any other value raises ValueError; the value has no effect with other strictness values.

complex_policy : str

What happens to an AND-linked subunit with no ortholog. Default "flag". "flag" keeps the reaction with a placeholder gene (OLD_<model>_<gene>) for later curator review, "keep" removes that subunit, "drop" drops the whole reaction. An isozyme group counts as a subunit without an ortholog only when none of its members has one, and then one placeholder stands for it. "drop" removes the reaction whenever a placeholder would remain, also when another isozyme branch is complete. Any other value raises ValueError.

only_genes_in_models : bool

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. As in RAVEN, it changes the mapping only with strictness 3: with 1 and 2 a template gene that no template contains is dropped in any case, and only the ValueError differs. Default False.

max_evalue : float

A hit passes only when its e-value is below this value. Default 1e-30, the RAVEN value. It makes no difference between 1e-4 and 1e-50, because the other two filters exclude every hit that it would exclude.

min_align_len : int

A hit passes only when its alignment length is at least this value. Default 100, as RAVEN's minLen. Values at or below 150 measured the same; on each of the four organisms tested, the score against KEGG annotations (wrong matches counted double) was 0.010 to 0.020 lower at 200 than at 100.

min_identity : float

A hit passes only when its identity is at least this value. Default 40. The value lies in the 35 to 45 range that the two reference sources favour; 35 scores higher for two organisms (yli 0.868 against 0.860, ani 0.814 against 0.804). It is the filter that limits the result.

review_identity : float or None

Lower identity threshold (for example 25) that collects the reactions the thresholds reject. Any reaction that transfers at this lower identity but not at min_identity is reported in HomologyResult.candidates with its strongest supporting hit and is not added to the model. A missing reaction can be gap-filled, whereas a wrongly transferred reaction is hard to find and remove, so the filters stay strict. It must be below min_identity, otherwise ValueError is raised. Hits with hits.attrs["source"] == "SuppOrth" have a placeholder identity of 100, so a warning is issued and no reaction is a candidate. Default None, which collects no candidates.

keep_gene_free : bool

Also transfers the template reactions that have no GPR (spontaneous, transport, exchange, artificial). By default they are dropped, as in RAVEN, because no gene supports them in the new organism. True transfers them, which suits a curated reference template whose gene-free reactions belong in every derived model. Default False.

preserve_notes : bool

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. Default False.

Returns

HomologyResult

The draft model, the template-to-new gene map and, when review_identity is given, the candidate reactions.

hits_from_supporth#

hits_from_supporth(source: str | Path | pd.DataFrame, organism_id: str, species_to_model: Mapping[str, str | None] | None=None, *, min_support: int=1, require_tools: Sequence[str]=(), min_identity: float | None=None, max_evalue: float | None=None, protein_to_gene: Mapping[str, str] | str | Path | None=None)

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 or 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 present in the mapping. 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: identity at least min_identity and e-value at most max_evalue. 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 transfer at the identity threshold review_identity but not at min_identity, with the best hit supporting each, sorted by identity, strongest first. None unless review_identity was given. These reactions are not in model.

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 fixed metrics (evalue 0, identity 100, align_len 1000, bitscore 1000, ppos 100). A known ortholog mapping can thus be passed to get_model_from_homology without a BLAST run.

Parameters

ortholog_pairs

Iterable of (source_gene, target_gene). The source gene belongs to the template model, the target gene to 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).

Parameters

path : str or Path

Path of the KEGG taxonomy file, plain or gzipped (suffix .gz).

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".

Parameters

path : str or Path

Path of the KEGG taxonomy file, plain or gzipped (suffix .gz).

domain : str

Domain name, or a prefix of it, for example "prokaryotes" or "euk". A category also matches when it is a prefix of domain. Raises ValueError when domain is empty.

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. Comment lines (#), empty lines and lines with fewer than 5 fields are skipped.

Parameters

text : str

Content of the --tblout file of hmmsearch.

parse_kegg_compounds#

parse_kegg_compounds(kegg_dir: str | Path)

Parse <kegg_dir>/compound and the optional compound.inchi into records.

Parameters

kegg_dir : str or Path

Directory of the arranged KEGG dump. It must contain the file compound and may contain compound.inchi. A compound with an InChI gets an empty formula.

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 or cobra YAML) plus the gzipped TSV tables into out_dir and returns {name: path} for everything written. The organism_gene_ko table is streamed to disk (see stream_organism_gene_ko) and not built in memory, so the full KEGG database can be processed; the other tables are built in memory.

When version is given (for example "kegg116"), the output filenames are version-prefixed (<version>_<name>), as in the published release assets; the keys of the returned dict are the logical table names.

Parameters

kegg_dir : str or Path

Directory of the arranged KEGG dump, as produced by download_kegg_dump.

out_dir : str or Path

Output directory. It is created, with parents, if missing.

version : str or None

Release tag, for example "kegg116", put before each file name as <version>_. When None (default), the file names have no prefix.

progress : bool

If True, report each stage and show a byte-progress bar over the organism_gene_ko streaming pass, the longest stage. Defaults to False.

parse_kegg_kos#

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

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

keep limits parsing to those KO ids (for example only KOs linked to reactions), as RAVEN's koList argument does; the gene lists are large.

Parameters

kegg_dir : str or Path

Directory of the arranged KEGG dump. It must contain the file ko.

keep : set of str or None

KO ids to parse. When None (default), all KOs are parsed.

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, set to irreversible for reactions that are irreversible across all KEGG maps (see _irreversible_from_mapformula).

Parameters

kegg_dir : str or Path

Directory of the arranged KEGG dump. It must contain the file reaction and may contain reaction_mapformula.lst.

parse_taxonomy#

parse_taxonomy(path: str | Path)

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

Parameters

path : str or Path

Path of the KEGG taxonomy file, plain or gzipped (suffix .gz).

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 (for example "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 and columns of phyl_dist.

Parameters

path : str or Path

Path of the KEGG taxonomy file, plain or gzipped (suffix .gz).

phyl_dist#

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

Compute the KEGG phylogenetic distance matrix from the taxonomy file.

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, or 1 if none coincide. This reproduces the values of RAVEN, including their asymmetry and negative values, so the output matches MATLAB keggPhylDist, which GECKO uses.

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

Time and memory use are O(N²), with N the number of KEGG organisms (~10⁴).

Parameters

path : str or Path

Path of the KEGG taxonomy file, plain or gzipped (suffix .gz).

only_in_kingdom : bool

If True, the distance between organisms of different domains is inf and the matrix has dtype float. If False (default), all distances are integers.

PhylDist#

KEGG phylogenetic distance structure of RAVEN getPhylDist and keggPhylDist.

ids and names are the KEGG organism ids and scientific names (names keep the trailing parenthetical of RAVEN; consumers that need them without it, for example geckopy, strip it). dist_matrix[i, j] is RAVEN's distance from ids[i] to ids[j]. It reproduces RAVEN's metric, which is asymmetric and may be negative; GECKO uses it only to select the closest (minimum-distance) organism, so the exact values matter for parity with MATLAB, and the sign and asymmetry do not affect that selection.

read_kegg_table#

read_kegg_table(path: str | Path)

Read a KEGG table from write_kegg_tables or stream_organism_gene_ko.

Compression is inferred from the file suffix. The published tables are gzipped TSV (.tsv.gz), including those with a version prefix (kegg116_<name>.tsv.gz). All columns are read as strings, and empty fields stay empty strings.

Parameters

path : str or Path

Path of the tab-separated table.

run_blast#

run_blast(organism_id: str, fasta: str | Path, model_ids: Sequence[str], ref_fastas: Sequence[str | Path], *, evalue: float=0.0001, threads: int | None=None, blastp: str | Path | None=None, makeblastdb: str | Path | None=None)

Bidirectional BLAST+ between an organism and template organisms.

Returns the hits DataFrame, filtered at evalue. Requires BLAST+ (blastp, makeblastdb). threads defaults to every CPU the process may use, as in RAVEN's getBlast.

The default evalue of 1e-4 equals the -evalue 10e-5 that MATLAB RAVEN's getBlast hardcodes. Here evalue is an argument. The cutoff is a pre-filter; the max_evalue, min_align_len and min_identity arguments of raven_toolbox.reconstruction.homology.homology.get_model_from_homology apply stricter limits downstream.

Parameters

organism_id : str

Identifier of the organism. It is written to the from_id column of the hits from the organism to a template and to the to_id column of the hits in the other direction.

fasta : str or Path

Protein FASTA file of the organism.

model_ids : Sequence[str]

Identifiers of the template models, written to the from_id and to_id columns. The sequence is parallel to ref_fastas and must have the same length, otherwise ValueError is raised.

ref_fastas : Sequence[str or Path]

Protein FASTA files of the templates. ref_fastas[i] is the reference proteome for template model_ids[i]. A single path is accepted as a list of one.

evalue : float

E-value cutoff passed to blastp. Default 1e-4.

threads : int or None

Number of threads passed to blastp. Default None, which uses every CPU the process may use.

blastp : str or Path or None

Path of the blastp binary. Default None, which resolves the binary with raven_toolbox.binaries.resolve_binary.

makeblastdb : str or Path or None

Path of the makeblastdb binary. Default None, which resolves the binary with raven_toolbox.binaries.resolve_binary.

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 | None=None, sensitivity: str='--more-sensitive', diamond: str | Path | None=None)

Bidirectional DIAMOND between an organism and template organisms.

Returns the hits DataFrame. Requires DIAMOND. threads defaults to every CPU the process may use, as in RAVEN's getDiamond.

The default evalue of 1e-3 equals DIAMOND's own default, which MATLAB RAVEN's getDiamond uses because it passes no --evalue. run_blast uses 1e-4, following getBlast. Neither cutoff has been calibrated against data, unlike the KO-assignment cutoffs of raven_toolbox.reconstruction.kegg.assign_kos, where measurement justified diverging from RAVEN.

Parameters

organism_id : str

Identifier of the organism. It is written to the from_id column of the hits from the organism to a template and to the to_id column of the hits in the other direction.

fasta : str or Path

Protein FASTA file of the organism.

model_ids : Sequence[str]

Identifiers of the template models, written to the from_id and to_id columns. The sequence is parallel to ref_fastas and must have the same length, otherwise ValueError is raised.

ref_fastas : Sequence[str or Path]

Protein FASTA files of the templates. ref_fastas[i] is the reference proteome for template model_ids[i]. A single path is accepted as a list of one.

evalue : float

E-value cutoff passed to diamond blastp. Default 1e-3.

threads : int or None

Number of threads passed to diamond blastp. Default None, which uses every CPU the process may use.

sensitivity : str

DIAMOND command-line flag appended unchanged to the diamond blastp command, for example "--more-sensitive" or "--sensitive". "" appends no flag and selects DIAMOND's default fast mode. Default "--more-sensitive".

diamond : str or Path or None

Path of the diamond binary. Default None, which resolves the binary with raven_toolbox.binaries.resolve_binary.

run_hmmsearch#

run_hmmsearch(fasta: str | Path, library: str | Path, *, threads: int | None=None, 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), which is the fast search direction. It runs with --cpu and needs no hmmpress. -Z is set to the profile count (when it is positive) so that the per-hit E-values match the convention assign_kos is calibrated against (identical to a hmmscan against the same library). A RuntimeError is raised if hmmsearch exits with a non-zero status.

Parameters

fasta : str or Path

Proteome FASTA file of the query organism (the search target).

library : str or Path

Concatenated multi-profile HMMER file of the KO library (the search query). It does not need hmmpress.

threads : int or None

Number of hmmsearch worker threads (--cpu). When None (default), every CPU the process may use, as in RAVEN's getKEGGModelForOrganism.

hmmsearch : str or Path or None

Path of the hmmsearch executable. When None (default), it is resolved by raven_toolbox.binaries.resolve_binary.

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 | None=None, extra_args: Sequence[str]=(), supporth: str | os.PathLike | None=None, **reader_kwargs)

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

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.

Requires the suppOrth command, found from supporth, then the RAVEN_SUPPORTH environment variable, then PATH. SuppOrth's predictors do not run on Windows. On Windows, SuppOrth can run under WSL or on a Linux host, and its suppOrth.tsv can be read with hits_from_supporth.

Parameters

organism_id : str

SuppOrth species label of the organism to build a model for. It is also the id of the organism in the returned hits. It must be non-empty, not "." or "..", have no leading or trailing space, contain none of ,, / and \, and differ from every id in model_ids; otherwise ValueError is raised.

fasta : str or Path

Protein FASTA file of the organism. A .gz file is read as gzip when the ids are checked.

model_ids : Sequence[str]

Ids of the template models, with the same restrictions as organism_id. The ids are passed to SuppOrth as --labels, so the table uses them as species labels. The sequence is parallel to ref_fastas and must have the same length, otherwise ValueError is raised. At least one template is needed.

ref_fastas : Sequence[str or Path]

Protein FASTA files of the templates. 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). A single path is accepted as a list of one.

outdir : str or Path

Output directory of SuppOrth. It is created when it does not exist. It receives suppOrth.tsv and run_supporth.log.

tools : str or Sequence[str]

Predictors to run, passed to SuppOrth as --tools: a comma-separated string or a sequence of names. "all" selects orthofinder, sonicparanoid and broccoli. Default "all".

support_lv : int or None

SuppOrth's --supportLv: a pair is kept if at least support_lv predictors found it, or if its identity is above identity_filter. It must be between 1 and 4, otherwise ValueError is raised. Default None, which uses 3, or the number of tools when fewer are selected. A warning is issued when it is larger than the number of tools and extra_args does not contain "--oma".

identity_filter : float

SuppOrth's --identity_filter in percent: a pair is kept if its identity is above this value, or if at least support_lv predictors found it. Default 80.0.

alignment : str

Aligner, passed to SuppOrth as --alignment. Default "diamond".

threads : int or None

Number of threads, passed to SuppOrth as -t. Default None, which uses every CPU the process may use, as in RAVEN's runSuppOrth.

extra_args : Sequence[str]

Arguments appended to the suppOrth run command, for example ["--resume"] or ["--oma", "path/to/Output"]. Default: none.

supporth : str or os.PathLike or None

Path of the suppOrth command. Default None, which uses the RAVEN_SUPPORTH environment variable, then suppOrth or supporth on PATH. A missing command raises FileNotFoundError.

**reader_kwargs : dict

Keyword arguments forwarded to hits_from_supporth (species_to_model, min_support, require_tools, min_identity, max_evalue and protein_to_gene).

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.

KEGG has ~9M gene and KO associations, which are 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 improves compression; the order also matches the by-organism query pattern in get_kegg_model_for_organism. The table is compressed with gzip, which MATLAB's built-in gunzip reads without an external tool; xz would give a smaller file.

The sort is an external merge sort: the row buffer is written as a sorted run, to a gzipped temporary file, once it reaches chunk_rows rows, and the runs are merged with heapq.merge, so peak memory does not grow with the size of KEGG. Only the ko_names table (one row per KO) is held in full; it is returned with the columns ko and name.

Parameters

kegg_dir : str or Path

Directory of the arranged KEGG dump. It must contain the file ko.

keep : set of str

KO ids to write, for example only the KOs linked to reactions, as keep of parse_kegg_kos.

ogk_path : str or Path

Path of the gzipped output table. Temporary sort files are written to its parent directory.

chunk_rows : int

Row count at which the buffer is written as a sorted run. It is checked after all GENES lines of a KO entry are added, so a run can exceed it by the rows of one entry. Defaults to 1,000,000.

progress : bool

If True, show a tqdm byte-progress bar over the (large) ko read pass. Defaults to False.

validate_hits#

validate_hits(hits: pd.DataFrame)

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

Parameters

hits : pandas.DataFrame

Hits table. It must contain every column in HIT_COLUMNS, otherwise ValueError is raised. Extra columns are allowed.

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 format shared with MATLAB RAVEN; MATLAB's built-in gunzip reads it without an external tool. prefix version-tags the filenames (for example kegg116_). Returns the written paths.

Parameters

tables : dict of str to pandas.DataFrame

Tables keyed by name, for example the result of build_kegg_tables. Each is written to <prefix><name>.tsv.gz without the index.

out_dir : str or Path

Output directory. It is created, with parents, if missing.

prefix : str

String put before each file name. Defaults to "".