Initialization (Python)#

raven-toolbox objects in raven_toolbox.init.

Functions#

Function

Summary

check_levels

Raise ValueError for a level in df that the upper-cased table lacks.

classify_reactions

Classify reactions into the ftINIT toIgnore categories (prepINITModel).

fill_tasks

Add minimum-cost reference reactions that make the tasks feasible in model.

ftinit

Run the full ftINIT pipeline on prep and return the context-specific model.

FtInitResult

Result of run_ftinit.

gene_scores_from_expression

Gene scores from RNA-seq or array expression, RAVEN's 5·ln(level/reference).

get_init_steps

Return the step schedule for a named ftINIT series (RAVEN getINITSteps).

group_rxn_scores

Sum original reaction scores over merged groups (RAVEN groupRxnScores).

hpa_gene_scores

Numeric gene scores from HPA levels for one tissue (optionally one celltype).

HPAData

Tidy HPA proteomics data: one row per (gene, tissue, cell type).

HPARnaData

Tidy HPA RNA-seq data: one row per (gene, tissue) with TPM.

InitStep

One ftINIT MILP step.

level_score_table

Level scores keyed by upper-cased name.

merge_linear

Merge linearly-dependent reactions; return (reduced, orig_ids, group_ids, reversed).

parse_hpa

Parse an HPA proteomics dump (normal_tissue.tsv; version 17 or later format).

parse_hpa_rna

Parse an HPA RNA-seq dump, as RAVEN's parseHPArna.

prep_init_model

Build PrepData from a template model, the work done once per template.

PrepData

One-time ftINIT preprocessing of a template model (RAVEN prepData).

ReactionMasks

Reaction-category id sets (RAVEN's toIgnore*), in 8-bit-pattern order.

remove_low_score_genes

Remove negative-scoring genes from GPRs (RAVEN removeLowScoreGenes).

rescale_for_init

Compress each reaction's stoichiometric dynamic range.

rna_gene_scores

Numeric gene scores from HPA RNA-seq TPM for one tissue.

run_ftinit

Run the single-step ftINIT MILP and return the extracted model.

score_for_init

Reaction scores and gene scores for tissue ({reaction id: score}, {gene id: score}).

score_reactions_from_genes

Score each reaction from per-gene scores via its GPR ({reaction_id: score}).

TaskFillResult

Result of fill_tasks: the gap-filled model and what was added.

TranscrData

Gene expression, as RAVEN's transcrData.

Reference#

check_levels#

check_levels(df: pd.DataFrame, table: Mapping[str, float])

Raise ValueError for a level in df that the upper-cased table lacks.

All rows count, whichever tissue is scored, as RAVEN's scoreModel checks every level of hpaData.levels.

classify_reactions#

classify_reactions(model: cobra.Model, *, ext_comp: str='e', spontaneous: Iterable[str]=(), custom: Iterable[str]=())

Classify reactions into the ftINIT toIgnore categories (prepINITModel).

A reaction is "GPR-less" when its gene rule is empty.

Parameters

model : cobra.Model

Model whose reactions are classified.

ext_comp : str

Keyword-only; default "e". Id of the extracellular compartment.

spontaneous : Iterable[str]

Keyword-only; default empty. Ids of the spontaneous reactions. Ids that are not reactions of model are ignored.

custom : Iterable[str]

Keyword-only; default empty. Ids of the reactions in the custom category. Ids that are not reactions of model are ignored.

fill_tasks#

fill_tasks(model: cobra.Model, reference_model: cobra.Model, tasks: Iterable[Task], *, rxn_scores: Mapping[str, float] | None=None, time_limit: float | None=_FILL_TIME_LIMIT, seed: int=_FILL_SEED, resolve_ties: bool=False, mutate_in_place: bool=False, verbose: bool=False)

Add minimum-cost reference reactions that make the tasks feasible in model.

Port of RAVEN fitTasks with gapFillMode 'preMerged' (the mode of ftINIT). Task by task, if a task is infeasible in the growing model, _gap_fill_task finds the minimum-cost set of reference reactions that restores it, and the function adds them; later tasks see the earlier additions. reference_model supplies the candidates (its reactions not yet in the model, excluding exchange and boundary reactions). Tasks with should_fail are ignored. Each gap-fill MILP is single-threaded, uses a fixed seed and is bounded by time_limit (RAVEN's 300 s).

Boundary reactions are closed while testing and solving each task, so the task inputs and outputs come solely from the ranged metabolite bounds of the task (RAVEN gap-fills the exchange-free model). The returned model keeps its boundary reactions. Tasks that could not be filled are listed in failed_tasks and reported with a warning; the model cannot perform those tasks.

Parameters

model : cobra.Model

Model to gap-fill. It is copied first unless mutate_in_place is True.

reference_model : cobra.Model

Model that supplies the candidate reactions. After each gap-fill, the Gurobi parameters Threads, Seed, IntFeasTol and TimeLimit (the generic solver timeout on other solvers) are restored. With resolve_ties=True on Gurobi, MIPGap and MIPGapAbs are restored as well.

tasks : Iterable[Task]

Tasks to make feasible in the result. Tasks with should_fail are skipped. A task that cannot be gap-filled stays infeasible, is listed in failed_tasks and gives a warning.

rxn_scores : Mapping[str, float] or None

Keyword-only; default None. Score of each reference reaction, keyed by original reaction id. The cost of a candidate is -min(score, -0.1), and a missing score gives cost 1. None gives every candidate cost 1.

time_limit : float or None

Keyword-only; default 300 s (RAVEN's value). Time limit of each gap-fill MILP in seconds. None sets no limit with Gurobi. At genome scale the tie-break phases of resolve_ties can exhaust it; the incumbent is then an unproven tie-break and a warning is issued, as in the main extraction.

seed : int

Keyword-only; default 26 (RAVEN's value). Random seed of each gap-fill MILP. It is applied to Gurobi only and ignored by other solvers.

resolve_ties : bool

Keyword-only; default False. When True, the degenerate min-cost fill is pinned to the fewest, lowest-id reactions, so that the added set does not depend on the solver seed or version (see _resolve_ties_fill).

mutate_in_place : bool

Keyword-only; default False. When True, model is gap-filled directly and returned as TaskFillResult.model; otherwise model is copied first and stays unchanged. Copying a genome-scale, Gurobi-backed model is expensive (see _gap_fill_task for the same effect on reference_model). raven_toolbox.init.ftinit passes True because its out is a disposable copy.

verbose : bool

Keyword-only; default False. Prints one line per task that is not skipped (added count and running total).

ftinit#

ftinit(prep, tissue: str | None=None, celltype: str | None=None, hpa: HPAData | None=None, *, transcr: TranscrData | HPARnaData | None=None, rxn_scores: Mapping[str, float] | None=None, gene_scores: Mapping[str, float] | None=None, remove_genes: bool=True, use_scores_for_tasks: bool=True, no_gene_score: float=-2.0, isozyme_scoring: str='max', complex_scoring: str='min', multiple_cell_scoring: str='max', hpa_level_scores: Mapping[str, float] | None=None, data_precedence: str='gene', series: str='1+1', steps=None, fill_gaps: bool=True, metabolomics: Iterable[str] | None=None, prod_weight: float=0.5, force_on: float=_FORCE_ON, big_m: float=_BIG_M, mip_gap: float | None=None, mip_gap_abs: float | None=10.0, time_limit: float | None=None, prove_abs_gap: float | None=None, resolve_ties: bool=False, seed: int=_EXTRACT_SEED, threads: int=_EXTRACT_THREADS, feas_tol: float=_EXTRACT_TOL, opt_tol: float=_EXTRACT_TOL, int_feas_tol: float=_EXTRACT_TOL, presolve: int=_EXTRACT_PRESOLVE, verbose: bool=False)

Run the full ftINIT pipeline on prep and return the context-specific model.

Like RAVEN's ftINIT(prepData, tissue, celltype, hpaData, 'transcrData', ...), the scores come from hpa, transcr or both for tissue (and celltype) on the reference model, by score_for_init. Alternatively pass rxn_scores and no data. Each step regroups the scores under its ignore_mask, fixes the reactions turned on by earlier steps as essential (in their flux direction) when the step uses 'essential', and solves run_ftinit on the merged model. Reactions never turned on (and not essential or left-in) are removed from the reference model; exchange reactions are always kept (RAVEN re-adds them).

If fill_gaps and prep carries tasks, reactions are added back to make the tasks feasible (raven_toolbox.init.fill_tasks); tasks with should_fail are skipped and a task that cannot be gap-filled stays infeasible, with a warning. When remove_genes is true and gene scores exist (given as gene_scores or derived from the data), negative-scoring genes are pruned from the GPRs at the end (raven_toolbox.init.remove_low_score_genes). The pruning aggregates with its own defaults ("max" for isozymes, "min" for complexes), whatever isozyme_scoring and complex_scoring are; RAVEN ftINIT calls removeLowScoreGenes without options too.

An essential reaction is forced to carry min(0.99·|carried flux|, force_on) in the forced direction, with force_on (default 0.1). The carried flux of every reaction starts at force_on and is replaced by the flux of each step after it is solved. This holds for the permanent (prep) essentials too. The per-reaction essential_force of run_ftinit sets the forced magnitude.

Parameters

prep : raven_toolbox.init.PrepData

The preprocessed template model.

tissue : str or None

Tissue to score for, matched case-insensitively. It is required when scoring from hpa or transcr and unused with rxn_scores. A ValueError is raised when it is None and rxn_scores is None.

celltype : str or None

Cell type to score for within tissue, matched case-insensitively. The default None does not restrict the cell type.

hpa : HPAData or None

Protein-level data. It is passed to score_for_init. A ValueError is raised when rxn_scores is given together with hpa or transcr.

transcr : TranscrData or HPARnaData or None

Expression data. It is passed to score_for_init.

rxn_scores : Mapping[str, float] or None

Original reaction id to score. When given, no scoring from data takes place, and gene_scores gives the gene scores for the final pruning.

gene_scores : Mapping[str, float] or None

Gene scores for the final pruning. They are used only when rxn_scores is given; otherwise the gene scores come from score_for_init.

remove_genes : bool

When True (default) and gene scores exist, low-score genes are pruned from the GPRs.

use_scores_for_tasks : bool

When True (default), the cost of the reactions added for the tasks is set from the reaction scores. When False, raven_toolbox.init.fill_tasks gets no scores.

no_gene_score : float

Score of a reaction without gene scores. The default is -2.0. It is passed to score_for_init.

isozyme_scoring : str

How the scores of OR alternatives in a GPR are combined: "min", "max", "median" or "average". The default is "max". It is passed to score_for_init.

complex_scoring : str

How the scores of AND subunits in a GPR are combined, with the same values as isozyme_scoring. The default is "min".

multiple_cell_scoring : str

How the data of several cell types or samples that match tissue are combined: "max" or "average". The default is "max".

hpa_level_scores : Mapping[str, float] or None

Maps an HPA expression level name to its score, compared case-insensitively. The default None uses RAVEN's hpaLevelScores.

data_precedence : str

"gene" (default) lets an HPA gene score replace an expression score of the same gene. "reaction" scores a reaction from HPA alone if any of its genes has an HPA score and from expression otherwise.

series : str

Names the staged schedule to run (see get_init_steps for what each option changes). The default is '1+1'. It is unused when steps is given.

steps : list of InitStep or None

A custom schedule, used in place of the one that series names. The default None uses get_init_steps(series).

fill_gaps : bool

When True (default) and prep carries tasks, reactions are added back to make the tasks feasible. Tasks with should_fail are skipped, and a task that cannot be gap-filled stays infeasible and gives a warning.

metabolomics : Iterable[str] or None

A list of detected metabolite names (case-insensitive, prep.ref_model metabolite names, matched and unioned across every metabolite sharing that name, for example the same compound in several compartments). Each is resolved once, before the step loop, into its producer reactions (translated through the merge groups, because the linear merge collapses them into combined reactions) and rewarded prod_weight in the objective of every step when any producer is on. The metabolomics entry of run_ftinit gives the MILP mechanics and its one caveat (the bonus checks that a producer is "on", not that the metabolite is net-produced through it). A name that matches no metabolite in prep.ref_model gives a warning. Ported from RAVEN 3's ftINITInternalAlg. The default None adds no metabolomics terms.

prod_weight : float

Reward per detected metabolite in the objective of every step. The default is 0.5 here, while RAVEN's ftINIT.m passes 5 to ftINITInternalAlg. It is unused when metabolomics is not given.

force_on : float

The flux that a reaction needs to count as "on", and the initial carried flux of every reaction. The default is 0.1.

big_m : float

Cap on the big-M of each scored reaction, which is otherwise its flux capacity. The default is 1000 (RAVEN's bigMCap). See run_ftinit and the module docstring.

mip_gap : float or None

Relative MIP gap for a step without a run schedule. A step with a schedule (every named series, including full) uses the gaps of its own runs (full: one run at 0.0004 and 10000 s). The default None leaves the solver default.

mip_gap_abs : float or None

Absolute gap, made relative to the objective, for a step without a run schedule; the default is 10.0. A step with a schedule uses the gaps of its own runs. On genome-scale models, the gaps and the time limit determine whether the MILP is tractable; the INIT parameter calibration study <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/init-param-calibration.md>_ on raven-docs gives the calibration table. Neither gap applies when prove_abs_gap is set.

time_limit : float or None

Caps the time limit of each run of every step, scheduled or not, in seconds. The default None applies the time limit of each scheduled run, or no limit for a step without a schedule. With prove_abs_gap set, the schedule is not read and the single solve of each step gets this value, so None sets no limit. The gap-fill of the tasks uses RAVEN's fixed per-task limit and ignores it.

prove_abs_gap : float or None

Opt-in; the default None gives RAVEN's behaviour. When set, replaces the loose relative-gap escalation with a single solve per step proven to this fixed absolute gap, so the kept-set objective is optimal at any scale and not an arbitrary within-gap incumbent. It changes which optimum is returned and not its score; the alternative optima are equally consistent with the expression data. The sweep ran on Human-GEM v2.0.0-38 with one cell line, four seeds (1234, 7, 1 and 2), resolve_ties=True, no gap-fill and time_limit 1800 s. The default escalation returned a step-1 objective 4.0 to 12.0 below the proven optimum; with seeds 1234 and 7 alone it ended 7.0 and 12.0 below, and 4.0 is seed 2. Every value from 2.0 to 0.05 reached the optimum for all four seeds and gave the same reaction sets, so 1.0 lies in the range of equal results; 5.0 stopped 2.0 below the optimum for one seed. Every step solve (the 48 primary solves) ended with status optimal; the flux phases phase2c-minflux (3 of 4 seeds) and phase2d-direction (4 of 4 seeds) reached the time limit and were adopted as unproven incumbents. The mean build time was about 5770 s at 0.05 against 5640 s at 1.0. The study recommends 1.0 together with resolve_ties when the seed-to-seed swing of resolve_ties alone (41 reactions, 1 gene here) is too large and the run time is acceptable. The full-pipeline runs of the calibration study used time_limit 600 s per step and the determinism builds 1800 s.

resolve_ties : bool

Opt-in; the default False gives RAVEN's behaviour. It reduces the arbitrariness of the tie-break in the degenerate MILP and gives a more parsimonious and more reproducible extracted model. It changes which optimum is returned and not its score. It adds a lexicographic phase 2 that selects the sparsest (then lowest-id) optimum, so the solver does not make the degenerate choice. It applies to each extraction step and to the task gap-fill. In each extraction step it then takes the fluxes from the least-total-flux solution for that on set, so the directions and forced magnitudes that the next step reads from them do not depend on the solver's vertex either; among those, it prefers the forward direction of the lower-ranked reaction, a merged reaction being ranked by the member with the smallest original id. The last step, whose fluxes no step reads, skips these flux phases. A step whose primary solve stops at time_limit keeps its primary solution without any of these phases.

seed : int

Gurobi's Seed (default 1234, RAVEN's value), forwarded to every step.

threads : int

Gurobi's Threads (default 1, RAVEN's value), forwarded to every step. More than one thread makes the result depend on thread timing.

feas_tol : float

Gurobi's FeasibilityTol (default 1e-9, RAVEN's value), forwarded to every step.

opt_tol : float

Gurobi's OptimalityTol (default 1e-9, RAVEN's value), forwarded to every step.

int_feas_tol : float

Gurobi's IntFeasTol (default 1e-9, RAVEN's value), forwarded to every step.

presolve : int

Gurobi's Presolve (default 2, RAVEN's value), forwarded to every step. Other values of the tolerances and of presolve change the optimum returned and the run time.

verbose : bool

When True, prints one line per staged step and is forwarded to fill_tasks for per-task gap-fill progress (as RAVEN ftINIT's console report). The default is False.

Notes

The studies on raven-docs ran at the default tolerances, presolve and one thread (see run_ftinit). The stability of the extraction under small input changes is a separate property that neither seed nor threads determines; see the ftINIT reproducibility study <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/ftinit-determinism.md>_ on raven-docs. Measured there on Human-GEM v2.0.0-38 with one cell line (seeds 1234 and 7, time_limit 1800 s): the default escalation differs between the seeds by 110 reactions and 17 predicted essential genes. resolve_ties alone gives identical reaction sets for repeated builds at a fixed seed and differs between the seeds by 41 reactions and 1 gene, at about 3.6 times the run time (about 14 min against 4 min). resolve_ties with prove_abs_gap=1.0 differs by 1 reaction and 0 genes, at about 26 times the baseline run time (about 100 min). The tie-break phases of resolve_ties can exhaust time_limit at genome scale; the spread is then reduced and not eliminated, and a warning is issued.

When models are compared before and after a curation, an edit applied to the template (followed by a new extraction) and the same edit applied to the extracted model are different controls: the new extraction contains the re-selection noise of the MILP in addition to the effect of the edit, while the edit applied to the extracted model changes only the reactions that it touches.

FtInitResult#

Result of run_ftinit.

gene_scores_from_expression#

gene_scores_from_expression(expression: Mapping[str, float], reference: Mapping[str, float] | float, *, factor: float=5.0, max_score: float=10.0, min_score: float=-5.0)

Gene scores from RNA-seq or array expression, RAVEN's 5·ln(level/reference).

reference is either a per-gene reference level (for example the cross-sample mean) or a single threshold for all genes: a gene expressed above its reference scores positive, below it negative. The score is factor * ln(level/reference), clamped to [min_score, max_score]. A zero level against a positive reference scores min_score; a positive level against a zero reference scores max_score; zero against zero, a negative level or reference, and a missing reference score min_score. A gene whose level is missing (None or NaN) gets no entry, so score_reactions_from_genes leaves it out of the rule.

Parameters

expression : Mapping[str, float]

Expression level of each gene, keyed by gene id.

reference : Mapping[str, float] or float

Per-gene reference level, or one threshold for all genes. A gene missing from a mapping has a missing reference.

factor : float

Keyword-only; default 5.0. Multiplies ln(level/reference).

max_score : float

Keyword-only; default 10.0. Upper limit of the score.

min_score : float

Keyword-only; default -5.0. Lower limit of the score.

get_init_steps#

get_init_steps(series: str='1+1', *, mets_to_ignore: Sequence[str]=())

Return the step schedule for a named ftINIT series (RAVEN getINITSteps).

'1+1' (default) has two steps; '2+1' has three steps, the first with positive reversibles dropped and metabolite secretion allowed; '1+0' and '2+0' omit the final step, which brings in the held-out GPR-less reactions; 'full' is a single MILP. mets_to_ignore are metabolite names removed from the stoichiometry in each step (for example H2O and H+); the 'full' series ignores none, as RAVEN's getINITSteps does. An unknown series raises ValueError.

Parameters

series : str

Name of the schedule, default "1+1". One of "1+1", "2+1", "1+0", "2+0" and "full"; any other value raises ValueError.

mets_to_ignore : Sequence[str]

Keyword-only; default empty. Names of the metabolites removed from the stoichiometry in each step. The "full" series ignores none, whatever the value.

group_rxn_scores#

group_rxn_scores(reduced_model: cobra.Model, orig_scores: Mapping[str, float], orig_rxn_ids: list[str], group_ids: list[int], to_zero: Iterable[str]=())

Sum original reaction scores over merged groups (RAVEN groupRxnScores).

orig_scores maps original reaction id to score (a missing id has score 0); to_zero are reactions to drop from the problem (the toIgnore masks), and their score becomes 0. A score of 0 and a group that cancels to 0 become 0.01, so that only ignore-masked reactions have a score of exactly 0 in the MILP. Returns {reduced_reaction_id: score}.

Parameters

reduced_model : cobra.Model

Reduced model from merge_linear. The result has one score per reaction of this model, and each of its reaction ids must be in orig_rxn_ids.

orig_scores : Mapping[str, float]

Score of each original reaction, keyed by reaction id. A missing id scores 0, and a score of 0 is set to 0.01 (see above).

orig_rxn_ids : list of str

Original reaction ids, in the order of group_ids.

group_ids : list of int

Group number of each original reaction, as returned by merge_linear; 0 marks a reaction that is not merged. It must have the length of orig_rxn_ids.

to_zero : Iterable[str]

Ids of original reactions whose score becomes 0 (the reactions dropped from the problem). The default is an empty tuple.

hpa_gene_scores#

hpa_gene_scores(hpa: HPAData, tissue: str, celltype: str | None=None, *, level_scores: Mapping[str, float] | None=None, multiple_celltype: str='best')

Numeric gene scores from HPA levels for one tissue (optionally one celltype).

Maps the categorical HPA levels to numbers via level_scores (default HPA_LEVEL_SCORES). As in RAVEN's scoreModel, level names, tissue and cell type are compared case-insensitively, and when the table has several names that differ only in case, the first one is used. Genes absent from the tissue are omitted from the output; score_reactions_from_genes then assigns no_gene_score to a reaction whose genes are all absent. A level that occurs anywhere in hpa and is not in the score table raises ValueError, whichever tissue is scored. A tissue or cell type that hpa does not contain gives an empty result. RAVEN raises only when neither the HPA nor the RNA data contain it; that check needs both data sets and is not made here.

Parameters

hpa : HPAData

Parsed HPA proteomics data, as returned by parse_hpa.

tissue : str

Tissue to score, compared case-insensitively with the tissue column of hpa.

celltype : str or None

Cell type to score, compared case-insensitively. The default None and an empty string use all cell types of the tissue.

level_scores : Mapping[str, float] or None

Keyword-only; default None, which uses HPA_LEVEL_SCORES; an empty mapping does too, as RAVEN's hpaLevelScores. Score of each HPA level name.

multiple_celltype : str

Keyword-only; default "best". Selects how the scores of several cell types of a tissue that carry the gene are combined: "best" is the maximum score (RAVEN default) and "average" is the mean over the cell types. Any other value raises ValueError.

HPAData#

Tidy HPA proteomics data: one row per (gene, tissue, cell type).

df columns: gene_id, gene_name, tissue, celltype, level, reliability. level is the categorical string from HPA; hpa_gene_scores maps it to numbers (a custom level_scores table can be passed to it).

HPARnaData#

Tidy HPA RNA-seq data: one row per (gene, tissue) with TPM.

df columns: gene_id, gene_name, tissue, tpm.

InitStep#

One ftINIT MILP step.

level_score_table#

level_score_table(level_scores: Mapping[str, float])

Level scores keyed by upper-cased name.

The first of the names that differ only in case is used.

merge_linear#

merge_linear(model: cobra.Model, no_merge: Iterable[str]=())

Merge linearly-dependent reactions; return (reduced, orig_ids, group_ids, reversed).

no_merge reaction ids are never merged. The reduced model carries no genes (merging makes GPRs meaningless); scores are remapped with group_rxn_scores.

Parameters

model : cobra.Model

Model to reduce. It is not modified.

no_merge : Iterable[str]

Ids of reactions that are never merged. The default is an empty tuple.

parse_hpa#

parse_hpa(path: str | Path)

Parse an HPA proteomics dump (normal_tissue.tsv; version 17 or later format).

The delimiter is detected automatically (HPA ships tab-separated files). The expected columns are Gene, Gene name, Tissue, Cell type, Level and Reliability; a missing column raises ValueError. Returns an HPAData with one row per (gene, tissue, cell type).

Parameters

path : str or Path

Path of the HPA proteomics file, read with pandas.read_csv.

parse_hpa_rna#

parse_hpa_rna(path: str | Path)

Parse an HPA RNA-seq dump, as RAVEN's parseHPArna.

One row per gene and tissue. The tissue column is Tissue (v19) or Sample (v18); the value column is TPM (v19) or Value (v18); Gene and Gene name are required. A file without a tissue or value column (for example one TPM column per tissue) raises ValueError. When a gene and tissue occur more than once, the last row is kept and a warning is issued. Rows whose value is not a number are dropped.

Parameters

path : str or Path

Path of the HPA RNA-seq file, read with pandas.read_csv.

prep_init_model#

prep_init_model(template: cobra.Model, tasks: Iterable[Task] | None=None, *, ext_comp: str='e', spontaneous: Iterable[str]=(), custom: Iterable[str]=(), essential_cache_path=None, simplify: bool=True, scale: bool=True, verbose: bool=False, processes: int=1)

Build PrepData from a template model, the work done once per template.

The result is shared by every ftINIT sample on this model.

With tasks, discovers the task-essential reactions (kept regardless of score), orients them irreversibly in their required direction, and drops tasks that are infeasible. Then classifies reactions into the omics-independent categories, linearly merges, and (unless scale=False) rescales the merged model's stoichiometry (rescale_for_init) so that a single MILP big-M is valid across all reactions.

Task-essential discovery runs with close_boundaries=True, as RAVEN's prepINITModel (which calls closeModel first), and the default of find_task_essential_reactions is additive. With close_boundaries=True the task file defines the whole boundary of the model during the discovery.

The two simplify_model passes are parallelised by cobra's own cobra.Configuration().processes.

The first simplification removes reactions only; every metabolite stays in the reference model while the tasks are checked, so a task can name any metabolite of the template, including ones whose names differ only in letter case. Afterwards the metabolites without a reaction are removed, except the ones that tasks name (inputs, outputs and equations, matched by exact name[compartment], failed tasks included; when several metabolites share a name[compartment], only the first in model order is kept). RAVEN prepINITModel does the same.

With the first simplification, every gene left without a reaction is removed from the reference model afterwards (RAVEN prepINITModel removes unused genes there), including genes that had no reaction in the template. Task ids may repeat (for example as category codes) and may be empty; a task that cannot be performed is removed by its position, so a task with the same id that can be performed stays. Raises ValueError when a task id is not text.

Parameters

template : cobra.Model

Template model. It is not modified; the preprocessing works on a copy.

tasks : Iterable[Task] or None

Metabolic tasks whose essential reactions are kept regardless of score. The default None skips the task-essential discovery.

ext_comp : str

Keyword-only; default "e". Id of the extracellular compartment, passed to classify_reactions.

spontaneous : Iterable[str]

Keyword-only; default empty. Ids of the spontaneous reactions, passed to classify_reactions.

custom : Iterable[str]

Keyword-only; default empty. Ids of the reactions in the custom category, passed to classify_reactions.

essential_cache_path : str or Path or None

Keyword-only; default None. Path of a file to which the accumulated results of the essential-reaction discovery are written, which makes the discovery resumable across interruptions; see find_task_essential_reactions. The default None writes no file.

simplify : bool

Keyword-only; default True. When False, the function skips RAVEN's dead-end and FVA-blocked removal and its reversibility-tightening passes, which leaves a larger merged model that may reach a different extraction optimum.

scale : bool

Keyword-only; default True. When False, the merged model's stoichiometry is not rescaled. In the INIT parameter calibration study on raven-docs (Human-GEM, one cell line), the full ftINIT pipeline without scaling was feasible at big_m 1000 and took 157 s against 83 s with scaling, and it was infeasible at big_m 100.

verbose : bool

Keyword-only; default False. Prints one line per phase plus per-task progress during essential-reaction discovery.

processes : int

Keyword-only; default 1. Number of processes for the task-essential discovery across tasks; see find_task_essential_reactions.

PrepData#

One-time ftINIT preprocessing of a template model (RAVEN prepData).

Built once per template, reused across samples. min_model is the merged model the MILP runs on; orig_rxn_ids/group_ids map its reactions back to the ref_model (the reference after the first simplification; the orientation of the essential reactions and the tightening of reversibility are undone on it after the merge, so its bounds and stoichiometry are the template's). essential_rxns are in merged ids and pre-oriented irreversibly (so the MILP forces flux forward). masks is on ref_model (= original) ids.

ReactionMasks#

Reaction-category id sets (RAVEN's toIgnore*), in 8-bit-pattern order.

ignored(pattern) returns the union of the categories whose bit is set, the reactions held out of that MILP step (they get score 0 and stay in the model).

remove_low_score_genes#

remove_low_score_genes(model: cobra.Model, gene_scores: Mapping[str, float], *, isozyme_scoring: str='max', complex_scoring: str='min', tied_reactions: list[str] | None=None)

Remove negative-scoring genes from GPRs (RAVEN removeLowScoreGenes).

Genes joined by OR are alternatives: the negative-scoring ones are removed, and when all alternatives are negative the least-negative one is kept. Genes joined by AND are subunits of a complex and are not removed individually.

Parameters

model : cobra.Model

The model to prune. It is not changed.

gene_scores : Mapping[str, float]

Gene id to score. Genes absent from it are unscored and never removed.

isozyme_scoring, complex_scoring : str

How the scores of OR alternatives and of complex subunits are aggregated ("max" and "min" by default).

tied_reactions : list of str, optional

When a list is given, the ids of the reactions in which an OR rule with only negative alternatives was cut to one of several alternatives with the same best score are appended to it, once per reaction. In such a rule the first of the tied alternatives in the rule is kept, so the result depends on the order of the rule. RAVEN keeps the first one too (CompMetEng/RAVEN#11).

Returns

tuple of (cobra.Model, list of str)

The pruned copy of the model, and the sorted ids of the genes that no reaction of the pruned copy uses; these genes are removed from the copy.

rescale_for_init#

rescale_for_init(model: cobra.Model, max_stoich_diff: float=25.0)

Compress each reaction's stoichiometric dynamic range.

A large spread of stoichiometric coefficients (for example a biomass or pool reaction with coefficients from 1e-3 to 1e2) forces correspondingly extreme flux magnitudes, so no single MILP big-M fits all reactions. As RAVEN, for each reaction the function caps every |coeff| at max_stoich_diff × min|coeff| (keeping signs), then scales the whole reaction so that its mean |coeff| is 1. Bounds are not changed. The function modifies model in place; only the merged MILP model is scaled, and the final output maps back to the unscaled ref_model by reaction id.

rna_gene_scores#

rna_gene_scores(rna: HPARnaData, tissue: str, *, reference: Mapping[str, float] | float | None=None, factor: float=5.0, max_score: float=10.0, min_score: float=-5.0)

Numeric gene scores from HPA RNA-seq TPM for one tissue.

Wrapper around raven_toolbox.init.score.gene_scores_from_expression (the same 5·ln(TPM/reference) scoring with clamping that is used elsewhere). The function selects the tissue, derives a reference if none is given (the per-gene mean TPM over all tissues, RAVEN's default for arrayData.threshold), and returns {gene_id: score}. As in RAVEN's scoreModel, the tissue is compared case-insensitively; when rna has several names that differ only in case, the highest TPM of a gene over those names is used (RAVEN's default multipleCellScoring='best'). A tissue that is not in rna raises ValueError.

Parameters

rna : HPARnaData

Parsed HPA RNA-seq data, as returned by parse_hpa_rna.

tissue : str

Tissue to score, compared case-insensitively with the tissue column of rna.

reference : Mapping[str, float] or float or None

Keyword-only. Per-gene reference TPM, or one threshold for all genes. The default None uses the mean TPM of each gene over all tissues.

factor : float

Keyword-only; default 5.0. Multiplies ln(TPM/reference).

max_score : float

Keyword-only; default 10.0. Upper limit of the gene score.

min_score : float

Keyword-only; default -5.0. Lower limit of the gene score.

run_ftinit#

run_ftinit(model: cobra.Model, rxn_scores: Mapping[str, float] | None=None, *, essential_rxns: Iterable[str] | None=None, essential_directions: Mapping[str, int] | None=None, essential_force: Mapping[str, float] | None=None, allow_excretion: bool=False, rem_pos_rev: bool=False, ignore_mets: Iterable[str]=(), force_on: float=_FORCE_ON, force_on_ess: float=_FORCE_ON, big_m: float=_BIG_M, mip_gap: float | None=None, mip_gap_abs: float | None=None, time_limit: float | None=None, prove_abs_gap: float | None=None, resolve_ties: bool=False, rank_keys: Mapping[str, str] | None=None, _pin_fluxes: bool=True, metabolomics: Mapping[str, Iterable[str]] | None=None, prod_weight: float=0.5, seed: int=_EXTRACT_SEED, threads: int=_EXTRACT_THREADS, feas_tol: float=_EXTRACT_TOL, opt_tol: float=_EXTRACT_TOL, int_feas_tol: float=_EXTRACT_TOL, presolve: int=_EXTRACT_PRESOLVE)

Run the single-step ftINIT MILP and return the extracted model.

The module docstring describes the formulation. raven_toolbox.init.ftinit calls this function once per step of its staged schedule.

Parameters

model : cobra.Model

The model to extract from. It is not changed.

rxn_scores : Mapping[str, float] or None

Reaction id to score. A reaction without a score has score 0 and is left free in the model, neither scored nor removable. The default None scores no reaction.

essential_rxns : Iterable[str] or None

Reaction ids forced to carry flux, by default none. The flux is at least force_on_ess, clamped to the reaction's capacity in the forced direction, that is ub forward and -lb reverse, and overridable per reaction via essential_force.

essential_directions : Mapping[str, int] or None

Maps an essential reaction id to +1 (forward) or -1 (reverse) for the forced direction. An essential reaction without an entry is forced forward.

essential_force : Mapping[str, float] or None

Maps an essential reaction id to the flux magnitude it is forced to carry, in place of force_on_ess. The default None uses force_on_ess for every essential reaction.

allow_excretion : bool

When True, the mass balance of each metabolite is relaxed from S·v = 0 to S·v >= 0. The default is False.

rem_pos_rev : bool

When True, a reversible reaction with a positive score gets score 0 and is handled as an unscored reaction, except that one which produces a metabolomics-detected metabolite still gets an indicator. The default is False.

ignore_mets : Iterable[str]

Metabolite names whose mass balance is dropped (RAVEN's per-step "simple metabolite" removal, for example H2O and H+). The default is empty.

force_on : float

Minimum absolute flux for a scored reaction to count as "on". The default is 0.1 (RAVEN's forceOnLim).

force_on_ess : float

Flux an essential reaction is forced to carry when essential_force has no entry for it. The default is 0.1.

big_m : float

Cap on the per-reaction big-M of a scored reaction. The default is 1000 (RAVEN's bigMCap). The module docstring describes how the big-M is derived.

mip_gap : float or None

Relative MIP gap. The default None leaves the solver default (Gurobi: MIPGap≈1e-4). It is applied on Gurobi only, and a named series of raven_toolbox.init.ftinit ignores it (every step has its own run schedule), so it acts on a step without a schedule. In a single step-0 solve of Human-GEM (one cell line; INIT parameter calibration study <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/init-param-calibration.md>_ on raven-docs) the solve took 24 to 25 s for every mip_gap from 0.0002 to 0.1, and the Jaccard index of the reaction set against the one from mip_gap 0.0002 fell below 0.97 from mip_gap 0.03.

mip_gap_abs : float or None

Absolute gap that is made relative to the objective. The default None disables it. When set and prove_abs_gap is None, an estimate solve is added before the final solve (limited to max(1, min(time_limit // 10, 30)) s when time_limit is set). On Gurobi, if the estimate objective obj is non-zero (step-2 objectives are negative), the final solve uses the relative gap min(max(mip_gap or 0, mip_gap_abs / |obj|), 1), which is RAVEN's AbsMIPGap/|objVal|.

time_limit : float or None

Time limit in seconds for each solve. The default None sets no limit and the solve runs until the gap is met.

prove_abs_gap : float or None

When set, proves the solve to this fixed absolute objective gap (on Gurobi: MIPGap 0 and MIPGapAbs this value); mip_gap and mip_gap_abs do not apply then. The default None leaves the gap to mip_gap and mip_gap_abs. In the four-seed sweep (seeds 1234, 7, 1 and 2; resolve_ties=True; no gap-fill; Human-GEM v2.0.0-38, one cell line) the default escalation returned a step-1 objective 4.0 to 12.0 below the proven optimum without a warning, while every value from 2.0 to 0.05 reached the optimum for all four seeds and gave the same reaction sets. Every step solve (the 48 primary solves) ended with status optimal; the flux phases phase2c-minflux (3 of 4 seeds) and phase2d-direction (4 of 4 seeds) reached the time limit and were adopted as unproven incumbents. The docstring of raven_toolbox.init.ftinit gives the measured run times.

resolve_ties : bool

Opt-in; the default False gives RAVEN's result. When True, the degeneracy of the MILP is resolved deterministically. After the score optimum is found, a lexicographic phase 2 holds the objective and minimises the number of kept negative-score reactions (phase 2a), then their summed id rank (phase 2b). The kept set is the sparsest optimum found, independent of the solver seed and version, where otherwise the solver returns an arbitrary tie-break. At genome scale the phase-2 solves can exhaust time_limit; the incumbent is then an unproven tie-break and a warning is issued. The tie-break is skipped when the primary solve itself stops at time_limit, because there is no proven optimum to hold. Phases 2a and 2b get half of time_limit each (rounded down, at least 1 s). After them, the reported fluxes are a least-total-flux solution for the selected on set (_min_abs_flux) and, among those, the one with the least id-rank-weighted backward flux (_min_backward_flux), each with the full time_limit. When one of these two flux phases finds no solution, the fluxes of the phase before it are kept and a warning is issued. See _resolve_ties.

rank_keys : Mapping[str, str] or None

Maps a reaction id to the key by which the id-rank phases of resolve_ties sort it (default: the id itself). raven_toolbox.init.ftinit passes the smallest original reaction id of each merge group, so the order does not depend on which member names the merged reaction.

metabolomics : Mapping[str, Iterable[str]] or None

In the id space of this model (raven_toolbox.init.ftinit resolves a list of metabolite names into this form once). It maps a detected-metabolite label to the set of reaction ids that produce it. The mon indicator of a metabolite (continuous, [0, 1]) is capped by the sum of its producers' own indicators and contributes prod_weight to the objective when any producer is on. This rewards keeping (or bringing back) a reaction that the data scored badly, because it makes a detected metabolite producible. A producer with score 0 is added to the problem (it gets a continuous indicator and contributes nothing itself) so that mon has an indicator to reference. The indicator of a negative-score producer gets a flux floor (as a positive reaction already has), so "on" means that the reaction carries flux. A detected metabolite is dropped without a warning if any of its producers is in essential_rxns, because it is produced regardless (RAVEN ftINITInternalAlg.m:106-111). The bonus requires that a producer is "on" and does not verify that the metabolite is net-produced through it: a reaction that consumes the metabolite on net can earn the bonus if it carries flux in reverse (a documented trade-off of RAVEN; the mass balance of the rest of the network makes the "on" flag meaningful). The default None adds no metabolomics terms.

prod_weight : float

Reward per metabolite in the mon terms of metabolomics. The default 0.5 is the default of RAVEN's ftINITInternalAlg when the argument is empty. It is unused when metabolomics is not given.

seed : int

Gurobi's Seed parameter (default 1234, RAVEN's value). The MILP is degenerate, so the seed selects which of many equal-score optima the solver returns. Varying it measures the dependence of a result on the tie-break. It is applied on Gurobi only.

threads : int

Gurobi's Threads (default 1, RAVEN's value). More threads can shorten the solve, and with more than one thread two runs of the same input can return different reaction sets, because multi-threaded Gurobi chooses between equal optima non-deterministically. The studies on raven-docs ran with one thread. The number of threads has no measured effect on the spread between different seeds, so it does not replace resolve_ties; see the ftINIT reproducibility study <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/ftinit-determinism.md>_ on raven-docs. It is applied on Gurobi only.

feas_tol : float

Gurobi's FeasibilityTol. The default is 1e-9 (RAVEN's value and Gurobi's minimum; Gurobi's own default is 1e-6). See the note on tolerances below.

opt_tol : float

Gurobi's OptimalityTol. The default is 1e-9 (Gurobi's own default is 1e-6).

int_feas_tol : float

Gurobi's IntFeasTol. The default is 1e-9 (Gurobi's own default is 1e-5). An "off" binary can have a value up to int_feas_tol and still permit flux up to its reaction's big-M times int_feas_tol, so a larger int_feas_tol loosens the on/off gate in proportion to that big-M.

presolve : int

Gurobi's Presolve (default 2, RAVEN's value).

Notes

feas_tol, opt_tol, int_feas_tol and presolve are RAVEN parity settings: they determine which optimal vertex the solver returns for the degenerate MILP and how binaries round at the 0.5 on/off cut. Changing them changes the optimum returned and the run time. The studies on raven-docs ran at the defaults; the effect of other values on the run time and on the extracted model has not been benchmarked.

The private parameter _pin_fluxes (default True) is set by raven_toolbox.init.ftinit: it is False in the last step, whose fluxes no later step reads. resolve_ties then skips phases 2c and 2d, and the reported fluxes are those of phase 2b.

score_for_init#

score_for_init(model: cobra.Model, hpa: HPAData | None, transcr: TranscrData | HPARnaData | None, tissue: str, celltype: str | None=None, *, no_gene_score: float=-2.0, isozyme_scoring: str='max', complex_scoring: str='min', multiple_cell_scoring: str='max', hpa_level_scores: Mapping[str, float] | None=None, data_precedence: str='gene')

Reaction scores and gene scores for tissue ({reaction id: score}, {gene id: score}).

Mirrors RAVEN's scoreModel. data_precedence="gene" lets an HPA score replace an expression score of the same gene; "reaction" scores a reaction from HPA alone if any of its genes has an HPA score and from expression otherwise. Genes without data are absent from the gene scores. Raises ValueError for missing data, unknown option values, a tissue or cell type found in neither dataset, or expression with fewer than two tissues and no threshold.

Parameters

model : cobra.Model

Model whose reactions are scored. Only genes of the model are scored.

hpa : HPAData or None

Parsed HPA proteomics data. None scores from transcr alone. At least one of hpa and transcr must be given.

transcr : TranscrData or HPARnaData or None

Expression data. An raven_toolbox.init.hpa.HPARnaData is converted with TranscrData.from_rna without a threshold. None scores from hpa alone.

tissue : str

Tissue to score, compared case-insensitively. It must occur in hpa or transcr.

celltype : str or None

Cell type to score, compared case-insensitively. The default None and an empty string use all cell types of the tissue. A non-empty value must occur in hpa or transcr.

no_gene_score : float

Keyword-only; default -2.0. Score of a reaction with no genes or no scored genes.

isozyme_scoring : str

Keyword-only; default "max". Combination of genes joined by OR; see raven_toolbox.init.score.score_reactions_from_genes.

complex_scoring : str

Keyword-only; default "min". Combination of genes joined by AND; accepts the same values as isozyme_scoring.

multiple_cell_scoring : str

Keyword-only; default "max". Combination of the values of several cell types (or samples) of the tissue for one gene: "max" takes the maximum and "average" the mean. Any other value raises ValueError.

hpa_level_scores : Mapping[str, float] or None

Keyword-only. Score of each HPA level name, compared case-insensitively. None and an empty mapping use RAVEN_HPA_LEVEL_SCORES.

data_precedence : str

Keyword-only; default "gene". "gene" lets an HPA gene score replace an expression score of the same gene. "reaction" scores a reaction from HPA alone if any of its genes has an HPA score and from expression otherwise. Any other value raises ValueError.

score_reactions_from_genes#

score_reactions_from_genes(model: cobra.Model, gene_scores: Mapping[str, float], *, isozyme_scoring: str='max', complex_scoring: str='min', no_gene_score: float=-2.0)

Score each reaction from per-gene scores via its GPR ({reaction_id: score}).

Isozymes (GPR OR) are combined with isozyme_scoring (default max); complex subunits (AND) with complex_scoring (default min). A reaction with no genes, or whose genes are all missing from gene_scores or scored NaN, gets no_gene_score (default -2).

Parameters

model : cobra.Model

Model whose reactions are scored.

gene_scores : Mapping[str, float]

Score of each gene, keyed by gene id.

isozyme_scoring : str

Keyword-only; default "max". Combination of genes joined by OR. One of "min", "max", "median" and "average"; any other value raises ValueError.

complex_scoring : str

Keyword-only; default "min". Combination of genes joined by AND. Accepts the same values as isozyme_scoring.

no_gene_score : float

Keyword-only; default -2.0. Score of a reaction with no genes or no scored genes.

TaskFillResult#

Result of fill_tasks: the gap-filled model and what was added.

TranscrData#

Gene expression, as RAVEN's transcrData.

levels is genes x samples (NaN where nothing was measured); tissues and celltypes name each column. threshold is a single value or one value per gene above which a gene counts as expressed; without it the mean over the measured samples is used.