Initialization (Python)#
raven-toolbox objects in raven_toolbox.init.
Functions#
Function |
Summary |
|---|---|
Raise |
|
Classify reactions into the ftINIT |
|
Add minimum-cost reference reactions that make the tasks feasible in |
|
Run the full ftINIT pipeline on |
|
Result of |
|
Gene scores from RNA-seq or array expression, RAVEN's |
|
Return the step schedule for a named ftINIT |
|
Sum original reaction scores over merged groups (RAVEN |
|
Numeric gene scores from HPA levels for one |
|
Tidy HPA proteomics data: one row per (gene, tissue, cell type). |
|
Tidy HPA RNA-seq data: one row per (gene, tissue) with TPM. |
|
One ftINIT MILP step. |
|
Level scores keyed by upper-cased name. |
|
Merge linearly-dependent reactions; return |
|
Parse an HPA proteomics dump ( |
|
Parse an HPA RNA-seq dump, as RAVEN's |
|
Build |
|
One-time ftINIT preprocessing of a template model (RAVEN |
|
Reaction-category id sets (RAVEN's |
|
Remove negative-scoring genes from GPRs (RAVEN |
|
Compress each reaction's stoichiometric dynamic range. |
|
Numeric gene scores from HPA RNA-seq TPM for one |
|
Run the single-step ftINIT MILP and return the extracted model. |
|
Reaction scores and gene scores for |
|
Score each reaction from per-gene scores via its GPR ( |
|
Result of |
|
Gene expression, as RAVEN's |
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.ModelModel whose reactions are classified.
ext_comp: strKeyword-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
modelare ignored.custom: Iterable[str]Keyword-only; default empty. Ids of the reactions in the custom category. Ids that are not reactions of
modelare 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.ModelModel to gap-fill. It is copied first unless
mutate_in_placeisTrue.reference_model: cobra.ModelModel that supplies the candidate reactions. After each gap-fill, the Gurobi parameters
Threads,Seed,IntFeasTolandTimeLimit(the generic solver timeout on other solvers) are restored. Withresolve_ties=Trueon Gurobi,MIPGapandMIPGapAbsare restored as well.tasks: Iterable[Task]Tasks to make feasible in the result. Tasks with
should_failare skipped. A task that cannot be gap-filled stays infeasible, is listed infailed_tasksand gives a warning.rxn_scores: Mapping[str, float] or NoneKeyword-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.Nonegives every candidate cost 1.time_limit: float or NoneKeyword-only; default 300 s (RAVEN's value). Time limit of each gap-fill MILP in seconds.
Nonesets no limit with Gurobi. At genome scale the tie-break phases ofresolve_tiescan exhaust it; the incumbent is then an unproven tie-break and a warning is issued, as in the main extraction.seed: intKeyword-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: boolKeyword-only; default
False. WhenTrue, 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: boolKeyword-only; default
False. WhenTrue,modelis gap-filled directly and returned asTaskFillResult.model; otherwisemodelis copied first and stays unchanged. Copying a genome-scale, Gurobi-backed model is expensive (see_gap_fill_taskfor the same effect onreference_model).raven_toolbox.init.ftinitpassesTruebecause itsoutis a disposable copy.verbose: boolKeyword-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.PrepDataThe preprocessed template model.
tissue: str or NoneTissue to score for, matched case-insensitively. It is required when scoring from
hpaortranscrand unused withrxn_scores. AValueErroris raised when it is None andrxn_scoresis None.celltype: str or NoneCell type to score for within
tissue, matched case-insensitively. The default None does not restrict the cell type.hpa: HPAData or NoneProtein-level data. It is passed to
score_for_init. AValueErroris raised whenrxn_scoresis given together withhpaortranscr.transcr: TranscrData or HPARnaData or NoneExpression data. It is passed to
score_for_init.rxn_scores: Mapping[str, float] or NoneOriginal reaction id to score. When given, no scoring from data takes place, and
gene_scoresgives the gene scores for the final pruning.gene_scores: Mapping[str, float] or NoneGene scores for the final pruning. They are used only when
rxn_scoresis given; otherwise the gene scores come fromscore_for_init.remove_genes: boolWhen True (default) and gene scores exist, low-score genes are pruned from the GPRs.
use_scores_for_tasks: boolWhen True (default), the cost of the reactions added for the tasks is set from the reaction scores. When False,
raven_toolbox.init.fill_tasksgets no scores.no_gene_score: floatScore of a reaction without gene scores. The default is -2.0. It is passed to
score_for_init.isozyme_scoring: strHow the scores of OR alternatives in a GPR are combined:
"min","max","median"or"average". The default is"max". It is passed toscore_for_init.complex_scoring: strHow the scores of AND subunits in a GPR are combined, with the same values as
isozyme_scoring. The default is"min".multiple_cell_scoring: strHow the data of several cell types or samples that match
tissueare combined:"max"or"average". The default is"max".hpa_level_scores: Mapping[str, float] or NoneMaps 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: strNames the staged schedule to run (see
get_init_stepsfor what each option changes). The default is'1+1'. It is unused whenstepsis given.steps: list of InitStep or NoneA custom schedule, used in place of the one that
seriesnames. The default None usesget_init_steps(series).fill_gaps: boolWhen True (default) and
prepcarries tasks, reactions are added back to make the tasks feasible. Tasks withshould_failare skipped, and a task that cannot be gap-filled stays infeasible and gives a warning.metabolomics: Iterable[str] or NoneA list of detected metabolite names (case-insensitive,
prep.ref_modelmetabolite 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 rewardedprod_weightin the objective of every step when any producer is on. Themetabolomicsentry ofrun_ftinitgives 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 inprep.ref_modelgives a warning. Ported from RAVEN 3'sftINITInternalAlg. The default None adds no metabolomics terms.prod_weight: floatReward per detected metabolite in the objective of every step. The default is 0.5 here, while RAVEN's
ftINIT.mpasses 5 toftINITInternalAlg. It is unused whenmetabolomicsis not given.force_on: floatThe flux that a reaction needs to count as "on", and the initial carried flux of every reaction. The default is 0.1.
big_m: floatCap on the big-M of each scored reaction, which is otherwise its flux capacity. The default is 1000 (RAVEN's
bigMCap). Seerun_ftinitand the module docstring.mip_gap: float or NoneRelative 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 NoneAbsolute 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 whenprove_abs_gapis set.time_limit: float or NoneCaps 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_gapset, 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 NoneOpt-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 andtime_limit1800 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 statusoptimal; the flux phasesphase2c-minflux(3 of 4 seeds) andphase2d-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 withresolve_tieswhen the seed-to-seed swing ofresolve_tiesalone (41 reactions, 1 gene here) is too large and the run time is acceptable. The full-pipeline runs of the calibration study usedtime_limit600 s per step and the determinism builds 1800 s.resolve_ties: boolOpt-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_limitkeeps its primary solution without any of these phases.seed: intGurobi's
Seed(default 1234, RAVEN's value), forwarded to every step.threads: intGurobi's
Threads(default 1, RAVEN's value), forwarded to every step. More than one thread makes the result depend on thread timing.feas_tol: floatGurobi's
FeasibilityTol(default 1e-9, RAVEN's value), forwarded to every step.opt_tol: floatGurobi's
OptimalityTol(default 1e-9, RAVEN's value), forwarded to every step.int_feas_tol: floatGurobi's
IntFeasTol(default 1e-9, RAVEN's value), forwarded to every step.presolve: intGurobi's
Presolve(default 2, RAVEN's value), forwarded to every step. Other values of the tolerances and ofpresolvechange the optimum returned and the run time.verbose: boolWhen True, prints one line per staged step and is forwarded to
fill_tasksfor per-task gap-fill progress (as RAVENftINIT'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 floatPer-gene reference level, or one threshold for all genes. A gene missing from a mapping has a missing reference.
factor: floatKeyword-only; default 5.0. Multiplies
ln(level/reference).max_score: floatKeyword-only; default 10.0. Upper limit of the score.
min_score: floatKeyword-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: strName of the schedule, default
"1+1". One of"1+1","2+1","1+0","2+0"and"full"; any other value raisesValueError.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.ModelReduced model from
merge_linear. The result has one score per reaction of this model, and each of its reaction ids must be inorig_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 strOriginal reaction ids, in the order of
group_ids.group_ids: list of intGroup number of each original reaction, as returned by
merge_linear; 0 marks a reaction that is not merged. It must have the length oforig_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: HPADataParsed HPA proteomics data, as returned by
parse_hpa.tissue: strTissue to score, compared case-insensitively with the
tissuecolumn ofhpa.celltype: str or NoneCell type to score, compared case-insensitively. The default
Noneand an empty string use all cell types of the tissue.level_scores: Mapping[str, float] or NoneKeyword-only; default
None, which usesHPA_LEVEL_SCORES; an empty mapping does too, as RAVEN'shpaLevelScores. Score of each HPA level name.multiple_celltype: strKeyword-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 raisesValueError.
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.ModelModel 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 PathPath 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 PathPath 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.ModelTemplate model. It is not modified; the preprocessing works on a copy.
tasks: Iterable[Task] or NoneMetabolic tasks whose essential reactions are kept regardless of score. The default
Noneskips the task-essential discovery.ext_comp: strKeyword-only; default
"e". Id of the extracellular compartment, passed toclassify_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 NoneKeyword-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; seefind_task_essential_reactions. The defaultNonewrites no file.simplify: boolKeyword-only; default
True. WhenFalse, 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: boolKeyword-only; default
True. WhenFalse, 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 atbig_m1000 and took 157 s against 83 s with scaling, and it was infeasible atbig_m100.verbose: boolKeyword-only; default
False. Prints one line per phase plus per-task progress during essential-reaction discovery.processes: intKeyword-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.ModelThe 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: strHow the scores of OR alternatives and of complex subunits are aggregated (
"max"and"min"by default).tied_reactions: list of str, optionalWhen 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: HPARnaDataParsed HPA RNA-seq data, as returned by
parse_hpa_rna.tissue: strTissue to score, compared case-insensitively with the
tissuecolumn ofrna.reference: Mapping[str, float] or float or NoneKeyword-only. Per-gene reference TPM, or one threshold for all genes. The default
Noneuses the mean TPM of each gene over all tissues.factor: floatKeyword-only; default 5.0. Multiplies
ln(TPM/reference).max_score: floatKeyword-only; default 10.0. Upper limit of the gene score.
min_score: floatKeyword-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.ModelThe model to extract from. It is not changed.
rxn_scores: Mapping[str, float] or NoneReaction 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 NoneReaction 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 isubforward and-lbreverse, and overridable per reaction viaessential_force.essential_directions: Mapping[str, int] or NoneMaps 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 NoneMaps an essential reaction id to the flux magnitude it is forced to carry, in place of
force_on_ess. The default None usesforce_on_essfor every essential reaction.allow_excretion: boolWhen True, the mass balance of each metabolite is relaxed from
S·v = 0toS·v >= 0. The default is False.rem_pos_rev: boolWhen 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: floatMinimum absolute flux for a scored reaction to count as "on". The default is 0.1 (RAVEN's
forceOnLim).force_on_ess: floatFlux an essential reaction is forced to carry when
essential_forcehas no entry for it. The default is 0.1.big_m: floatCap 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 NoneRelative 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.ftinitignores 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 everymip_gapfrom 0.0002 to 0.1, and the Jaccard index of the reaction set against the one frommip_gap0.0002 fell below 0.97 frommip_gap0.03.mip_gap_abs: float or NoneAbsolute gap that is made relative to the objective. The default None disables it. When set and
prove_abs_gapis None, an estimate solve is added before the final solve (limited tomax(1, min(time_limit // 10, 30))s whentime_limitis set). On Gurobi, if the estimate objectiveobjis non-zero (step-2 objectives are negative), the final solve uses the relative gapmin(max(mip_gap or 0, mip_gap_abs / |obj|), 1), which is RAVEN's AbsMIPGap/|objVal|.time_limit: float or NoneTime 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 NoneWhen set, proves the solve to this fixed absolute objective gap (on Gurobi:
MIPGap0 andMIPGapAbsthis value);mip_gapandmip_gap_absdo not apply then. The default None leaves the gap tomip_gapandmip_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 statusoptimal; the flux phasesphase2c-minflux(3 of 4 seeds) andphase2d-direction(4 of 4 seeds) reached the time limit and were adopted as unproven incumbents. The docstring ofraven_toolbox.init.ftinitgives the measured run times.resolve_ties: boolOpt-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 attime_limit, because there is no proven optimum to hold. Phases 2a and 2b get half oftime_limiteach (rounded down, at least 1 s). After them, the reportedfluxesare 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 fulltime_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 NoneMaps a reaction id to the key by which the id-rank phases of
resolve_tiessort it (default: the id itself).raven_toolbox.init.ftinitpasses 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 NoneIn the id space of this model (
raven_toolbox.init.ftinitresolves a list of metabolite names into this form once). It maps a detected-metabolite label to the set of reaction ids that produce it. Themonindicator of a metabolite (continuous,[0, 1]) is capped by the sum of its producers' own indicators and contributesprod_weightto 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 thatmonhas 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 inessential_rxns, because it is produced regardless (RAVENftINITInternalAlg.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: floatReward per metabolite in the
monterms ofmetabolomics. The default 0.5 is the default of RAVEN'sftINITInternalAlgwhen the argument is empty. It is unused whenmetabolomicsis not given.seed: intGurobi's
Seedparameter (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: intGurobi'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 replaceresolve_ties; see theftINIT 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: floatGurobi'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: floatGurobi's
OptimalityTol. The default is 1e-9 (Gurobi's own default is 1e-6).int_feas_tol: floatGurobi's
IntFeasTol. The default is 1e-9 (Gurobi's own default is 1e-5). An "off" binary can have a value up toint_feas_toland still permit flux up to its reaction's big-M timesint_feas_tol, so a largerint_feas_tolloosens the on/off gate in proportion to that big-M.presolve: intGurobi'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.ModelModel whose reactions are scored. Only genes of the model are scored.
hpa: HPAData or NoneParsed HPA proteomics data.
Nonescores fromtranscralone. At least one ofhpaandtranscrmust be given.transcr: TranscrData or HPARnaData or NoneExpression data. An
raven_toolbox.init.hpa.HPARnaDatais converted withTranscrData.from_rnawithout a threshold.Nonescores fromhpaalone.tissue: strTissue to score, compared case-insensitively. It must occur in
hpaortranscr.celltype: str or NoneCell type to score, compared case-insensitively. The default
Noneand an empty string use all cell types of the tissue. A non-empty value must occur inhpaortranscr.no_gene_score: floatKeyword-only; default -2.0. Score of a reaction with no genes or no scored genes.
isozyme_scoring: strKeyword-only; default
"max". Combination of genes joined by OR; seeraven_toolbox.init.score.score_reactions_from_genes.complex_scoring: strKeyword-only; default
"min". Combination of genes joined by AND; accepts the same values asisozyme_scoring.multiple_cell_scoring: strKeyword-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 raisesValueError.hpa_level_scores: Mapping[str, float] or NoneKeyword-only. Score of each HPA level name, compared case-insensitively.
Noneand an empty mapping useRAVEN_HPA_LEVEL_SCORES.data_precedence: strKeyword-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 raisesValueError.
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.ModelModel whose reactions are scored.
gene_scores: Mapping[str, float]Score of each gene, keyed by gene id.
isozyme_scoring: strKeyword-only; default
"max". Combination of genes joined by OR. One of"min","max","median"and"average"; any other value raisesValueError.complex_scoring: strKeyword-only; default
"min". Combination of genes joined by AND. Accepts the same values asisozyme_scoring.no_gene_score: floatKeyword-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.