Analysis (Python)

raven-toolbox objects in raven_toolbox.analysis.

Functions

Function	Summary
analyze_sampling	Compare the significance of a flux change with a gene-expression change.
compare_fluxes	Find the reactions whose flux differs between two conditions.
CompareFluxesResult	Outcome of a compare_fluxes comparison.
find_good_reactions	Reactions usable as random objectives: those that carry non-loop flux.
FluxSamplingResult	Output of raven_toolbox.analysis.sampling.random_sampling.
FluxWalker	Steppable navigator over a solved flux distribution.
fseof	Run FSEOF for over-production of target_rxn's product.
FSEOFResult	FSEOF output.
get_flux_z	Z-score for how much flux changed between two random-sampling runs.
get_min_nr_fluxes	Find a flux distribution activating as few reactions as possible.
max_volume_ellipsoid	Maximum-volume ellipsoid inscribed in the polytope {z : A z <= b}.
MetaboliteGroup	One metabolite the current reaction touches, and its other carriers.
MinNrFluxesResult	Outcome of a minimum-cardinality flux search.
NeighborReaction	One neighbouring reaction, as displayed under a metabolite group.
print_flux_path	Print trace_flux_path's result and return it.
random_sampling	Sample model's flux space with any of the sampling methods.
reporter_metabolites	Compute Reporter Metabolites from per-gene differential-expression p-values.
ReporterResult	Reporter-metabolite scores for one gene set.
trace_flux_path	Find the highest-flux-fraction path from from_rxn to to_rxn.
TraceFluxPathResult	Outcome of a trace_flux_path search.
walk_fluxes	Interactively navigate a flux distribution reaction by reaction.

Reference

analyze_sampling

analyze_sampling(t_expression: Sequence[float] | np.ndarray | pd.Series, df: float, solutions_a: pd.DataFrame, solutions_b: pd.DataFrame)

Compare the significance of a flux change with a gene-expression change.

Parameters

t_expression:

t-score for the change in gene expression, one per reaction, in the same order as solutions_a's columns. Positive means expression increased, negative means it decreased. A p-value can be converted to this by the inverse error function, with the sign chosen by hand because the p-value does not carry direction.

df:

Degrees of freedom for the expression t-test.

solutions_a, solutions_b:

Sampled flux vectors for the reference and test condition, samples × reactions; see get_flux_z.

Returns

pd.DataFrame

Indexed by reaction id, with columns:

• concordant: probability that flux and expression both changed significantly, in the same direction (RAVEN's pR).

• expression_only: probability that expression changed significantly but flux did not, or that they changed in opposite directions (pH).

• flux_only_or_discordant: probability that flux changed significantly but expression did not, or that they changed in opposite directions (pM).

Notes

Two differences from analyzeSampling.m. First, analyzeSampling.m computes an "opposite flux direction" branch (pM=erf(|Z|), pH=(1-pM)*..., pR=0) whose results are unconditionally overwritten by the next branch before scores is read; this dead code is not replicated. Second, analyzeSampling.m bundles its own tcdf (a pure-MATLAB fallback for when the Statistics Toolbox is not installed); this function uses scipy.stats.t.cdf.

compare_fluxes

compare_fluxes(model: cobra.Model, reference: pd.Series | Mapping[str, float], comparison: pd.Series | Mapping[str, float], *, cutoff: float=1e-08, metabolite_list: Sequence[str] | None=None)

Find the reactions whose flux differs between two conditions.

Parameters

model:

The model both flux vectors are for.

reference:

Flux for the reference condition (for example wild type), indexed by reaction id. Every reaction in the model must be present.

comparison:

Flux for the condition being compared against it (for example a mutant).

cutoff:

A reaction is active when |flux| > cutoff, and changed when |flux2 - flux1| > cutoff. Both comparisons are strict.

metabolite_list:

Metabolite names, compared case-insensitively against Metabolite.name. If given, only reactions touching at least one of them are compared; every other reaction is omitted from the result. A name matching no metabolite is reported in CompareFluxesResult.missing_metabolites and does not raise.

Returns

CompareFluxesResult :

CompareFluxesResult

Outcome of a compare_fluxes comparison.

changed holds one row per reaction whose |flux2 - flux1| exceeded the cutoff, sorted by abs_delta descending, with columns reaction, name, flux1, flux2, abs_delta, rel_change and type. rel_change is (flux2 - flux1) / |flux1|, and is NaN when |flux1| is below the cutoff, because the reaction has no baseline. type is "on", "off", "flip" or "".

turned_on, turned_off and flipped are reaction-id lists in the model's own order, not the sorted order of changed. missing_metabolites lists the names of metabolite_list that match no metabolite, and n_considered is the number of reactions compared.

find_good_reactions

find_good_reactions(model: cobra.Model, *, flux_tol: float=1e-09, loopless: bool=True, exclude_reactions: Iterable[str] | None=None)

Reactions usable as random objectives: those that carry non-loop flux.

A reaction is kept if the larger of |minimum| and |maximum| of its FVA range exceeds flux_tol. With loopless (default) the FVA is loopless (cycleFreeFlux), so reactions that can carry flux only through a thermodynamically infeasible cycle have a ~0 loopless range and are dropped. A fixed bound threshold would also drop legitimate reactions that reach the model's default (for example 1000) bound. With loopless=False the FVA is faster and keeps any flux-carrying reaction (loops included).

Parameters

model : cobra.Model

Model whose reactions are tested.

flux_tol : float

A reaction is kept when the larger of |minimum| and |maximum| of its FVA range exceeds this value. The default is 1e-9.

loopless : bool

Whether the FVA is loopless (cycleFreeFlux). The default is True.

exclude_reactions : Iterable[str] or None

Identifiers of reactions dropped from the result. The default None excludes none.

FluxSamplingResult

Output of raven_toolbox.analysis.sampling.random_sampling.

Attributes

samples:

Flux vectors shaped n_samples × n_reactions (one sample per row, reaction ids as columns; the cobra.sampling layout).

method:

"achr", "chrr", or "random_objective".

good_reactions:

random_objective only: reaction ids eligible as random objectives; None for the MCMC methods.

n_dimensions:

chrr only: dimension of the full-dimensional flux polytope sampled (degrees of freedom after fixing implicitly-determined reactions).

mve_converged:

chrr only: whether the MVE rounding solver reached its tolerance. With False, the last ellipsoid iterate is used as the rounding; elongated polytopes may need more thinning.

n_warmup:

achr only: number of FVA warmup directions.

fixed_reactions:

chrr only: reactions added to the equality system as fixed (FVA width below fixed_width_tol).

FluxWalker

Steppable navigator over a solved flux distribution.

fluxes is indexed by reaction id (for example a cobra.Solution's .fluxes, or the result of pfba(model).fluxes). Only a metabolite whose net contribution from the current reaction is at least cutoff in absolute value is shown, and at most max_per_met of its other flux-carrying reactions (> cutoff, ranked by their own |flux|) are listed per metabolite.

fseof

fseof(model: cobra.Model, target_rxn: str, *, biomass_rxn: str | None=None, n_steps: int=10, max_fraction: float=0.9, correlation_threshold: float=0.9, flux_eps: float=1e-06, min_target: float | None=None, equality: bool=False)

Run FSEOF for over-production of target_rxn's product.

Enforces the target flux in n_steps evenly spaced levels, from max_fraction / n_steps times the theoretical maximum up to max_fraction times the theoretical maximum (the target flux of the pFBA solution that maximises it), maximising growth (biomass_rxn or the model's current objective) with pFBA at each step. Returns an FSEOFResult.

Parameters

model : cobra.Model

Model to scan. It is restored after the call.

target_rxn : str

Identifier of the reaction whose flux is enforced. The theoretical maximum is found by maximising this reaction; a ValueError is raised when that maximum is not above flux_eps.

biomass_rxn : str or None

Identifier of the reaction maximised at each enforced level. The default None keeps the model's current objective.

n_steps : int

Number of enforced flux levels. The default is 10. Levels are evenly spaced; the scan raises RuntimeError when fewer than two levels are feasible.

max_fraction : float

Fraction of the theoretical maximum target flux at which the highest level is set. The default is 0.9.

correlation_threshold:

Minimum absolute Pearson correlation between the absolute flux of a reaction and the enforced product flux for the reaction to be reported as a target; a higher value keeps fewer reactions, those with the most consistent trend.

flux_eps:

Tolerance below which the theoretical maximum of the target flux is treated as zero (the scan then raises). It gates only that check.

min_target:

Lowest enforced flux. The n_steps levels are evenly spaced from min_target up to the scan's ceiling (max_fraction times the theoretical maximum). Defaults to None, in which case the lowest level is the ceiling divided by n_steps. Passing the flux of the target reaction at the growth optimum starts the scan at that flux.

equality:

How the enforced flux is imposed. False (the default) sets it as the lower bound of the target reaction, as RAVEN's FSEOF does; True fixes the target flux at the level (lower and upper bound), as GECKO's ecFSEOF does.

Notes

Reaction fluxes are classified as in RAVEN's FSEOF: a reaction is a target when some absolute flux exceeds 1e-10, the standard deviation of its absolute flux is above 0 and the correlation reaches correlation_threshold; it is knockout when the slope is negative and the last absolute flux is at most 1e-8, knockdown for a negative slope otherwise.

Each level is solved in two steps as in RAVEN: the objective is maximised, then the flux sum is minimised with the objective held at 0.999999 of its optimum at that level. A level at which the first step is infeasible ends the scan.

The top of the scan is max_fraction of the target flux in the pFBA solution that maximises the target, as in RAVEN's FSEOF and GECKO's ecFSEOF.

Differences from FSEOF.m: the scan stops at the first infeasible enforced level, so the last absolute flux used for the knockout label is that of the last feasible level; targets is sorted by |slope| * |r| (equal scores in model order) and is one table, whereas FSEOF.m prints the targets grouped by the labels Amplify, Knockdown and Knockout; it is not split into reaction and transport targets.

FSEOFResult

FSEOF output.

scan is reactions × enforced-flux-levels (the full flux scan); enforced are the enforced target fluxes; targets is the classified per-reaction table (sorted by score). gene_targets aggregates targets to genes.

get_flux_z

get_flux_z(solutions_a: pd.DataFrame, solutions_b: pd.DataFrame)

Z-score for how much flux changed between two random-sampling runs.

Parameters

solutions_a:

Sampled flux vectors for the reference condition, samples × reactions (for example random_sampling(...).samples).

solutions_b:

Sampled flux vectors for the test condition, same reaction columns.

Returns

pd.Series

Z-score per reaction (indexed by solutions_a's columns): positive when flux increased from A to B, negative when it decreased, 0 when the means are equal, clamped to [-100, 100]. A reaction with zero variance in both conditions and different means gets ±100.

Notes

A flux decrease with no variance in either condition scores -100, the same sign the general formula gives, as in RAVEN's getFluxZ.

get_min_nr_fluxes

get_min_nr_fluxes(model: cobra.Model, to_minimize: Sequence[str] | None=None, *, scores: Sequence[float] | None=None, big_m: float | None=None)

Find a flux distribution activating as few reactions as possible.

Parameters

model:

Model to solve on, subject to its own steady-state constraints and reaction bounds as given (the model's objective plays no role; a minimal solution among optimal ones requires constraining the objective beforehand, for example pinning it near its optimum).

to_minimize:

Reaction ids to give a binary "is this active" indicator. Defaults to every reaction in the model.

scores:

Per-reaction weights, same length and order as to_minimize. Negative scores discourage that reaction from being active (the more negative, the stronger); entries of 0 or above are replaced by the largest (least negative) negative score present, or by 0 if scores has no negative entries, as RAVEN's getMinNrFluxes does. Defaults to an implicit -1 for every reaction (unweighted cardinality minimisation).

big_m:

Override the big-M constant. Defaults to an estimate from the model's own bounds (see module docstring).

Returns

MinNrFluxesResult :

max_volume_ellipsoid

max_volume_ellipsoid(A: np.ndarray, b: np.ndarray, x0: np.ndarray | None=None, *, maxiter: int=150, tol: float=1e-06, reg: float=1e-08)

Maximum-volume ellipsoid inscribed in the polytope {z : A z <= b}.

Solves max log det E over centre x and SPD E such that the ellipsoid {x + E s : ||s||_2 <= 1} is contained in {z : A z <= b}, using the Zhang & Gao (2003) primal-dual interior-point method (the regularised variant shipped in COBRA's chrrSampler).

Parameters

A, b:

Polytope A z <= b with A of shape (m, n), m >= n + 1, and a non-empty interior.

x0:

A strictly interior point (A x0 < b). If None a Chebyshev centre is computed by LP.

maxiter, tol, reg:

Iteration cap, convergence tolerance, and the diagonal (Levenberg) regularisation added to the two Newton solves.

Returns

x:

Ellipsoid centre, shape (n,).

E:

Lower-triangular rounding transform with E @ E.T == E2 (the SPD ellipsoid matrix). The rounding substitution is z = x + E y.

converged:

Whether tol was reached within maxiter.

Notes

Validated against analytic cases: a box maps to the unit ball (E = I); a scaled or sheared box {-1 <= M z <= 1} gives E E^T = M^{-1} M^{-T}; the standard simplex gives its inscribed ball.

MetaboliteGroup

One metabolite the current reaction touches, and its other carriers.

MinNrFluxesResult

Outcome of a minimum-cardinality flux search.

Parameters

fluxes:

Flux for every reaction in the model (empty if infeasible).

active:

Ids from to_minimize whose indicator was on in the solution.

status:

Solver status string (for example "optimal", "infeasible").

NeighborReaction

One neighbouring reaction, as displayed under a metabolite group.

print_flux_path

print_flux_path(model: cobra.Model, fluxes: pd.Series | Mapping[str, float], from_rxn: str, to_rxn: str, *, print_fn=print, **kwargs)

Print trace_flux_path's result and return it.

The output is a simplified rendering of traceFluxPath.m's console diagram. Parity with RAVEN is defined on the structured result.

Parameters

model : cobra.Model

Model both reactions belong to.

fluxes : pd.Series or Mapping[str, float]

Flux by reaction id, as for trace_flux_path.

from_rxn : str

Source reaction id.

to_rxn : str

Target reaction id.

print_fn : callable

Function called with each output line. The default is the builtin print.

**kwargs : dict

Keyword arguments forwarded to trace_flux_path: cutoff, max_hops, trace_material, carbon_only and exclude_mets.

random_sampling

random_sampling(model: cobra.Model, n_samples: int=1000, *, method: str='achr', seed: int | None=None, thinning: int=100, warmup: int=1000, fixed_width_tol: float=1e-07, n_objectives: int=2, good_reactions: Iterable[str] | None=None, replace_max_bound: bool=False, min_flux: bool=False, loopless_good_reactions: bool=True, exclude_reactions: Iterable[str] | None=None, max_attempts: int=100, suppress_errors: bool=False, n_proc: int | None=None)

Sample model's flux space with any of the sampling methods.

Dispatches on method:

• "achr" (default): Artificially Centered Hit-and-Run; near-uniform MCMC sampling of the polytope interior (wraps cobra.sampling.ACHRSampler).

• "chrr": Coordinate Hit-and-Run with Rounding; near-uniform MCMC with maximum-volume-ellipsoid rounding, for better mixing on thin or ill-conditioned polytopes such as enzyme-constrained models.

• "random_objective": the random-objective vertex method of Bordel et al. (2010): each sample maximises a small random objective, returning a polytope vertex.

The "achr" and "chrr" methods draw the (near-)uniform interior distribution; "random_objective" draws diverse vertices. The samples are conditioned on the constraints set on the model before the call (for example a biomass lower bound, measured fluxes, enzyme-usage bounds).

Parameters

model : cobra.Model

Model to sample. Its current bounds and constraints define the sampled space.

n_samples:

Number of flux vectors to return.

method:

"achr" (default), "chrr", or "random_objective".

seed:

Seed for reproducible chains and draws.

thinning:

achr and chrr: Markov-chain steps between recorded samples (default 100). The default is calibrated for small models (~200 to 500 reactions). On genome-scale models (>2000 reactions) lag-1 autocorrelation remains high even at thinning=100 (measured: 0.926 on yeast-GEM with 4102 reactions), giving about 12 effective samples from 300 stored samples (ESS ≈ n × (1−ρ)/(1+ρ)). Independent chains at these default settings also fail to agree with each other: on yeast-GEM, 67.5% of reactions fail the Gelman-Rubin R-hat>1.1 convergence threshold across 4 independent chains (96.5% fail the stricter R-hat>1.01). Default-settings genome-scale output is therefore unconverged for most reactions. Options for genome-scale analyses are increasing thinning substantially (≥1000), increasing n_samples to compensate, using method='chrr' (raven-toolbox's rounding-based alternative to ACHR), or assessing sample quality with an ESS or R-hat diagnostic (see the sampling convergence calibration study <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/sampling-convergence-calibration.md>_ on raven-docs).

warmup:

chrr: burn-in steps discarded before the first recorded sample.

fixed_width_tol:

chrr: a reaction whose FVA range is narrower than this is added to the equality system as fixed, which keeps the reduced polytope full-dimensional.

n_objectives, good_reactions, loopless_good_reactions, exclude_reactions, max_attempts, suppress_errors:

random_objective only: see the method's parameters; good_reactions can be passed back from a previous result to skip the one-off FVA.

replace_max_bound:

random_objective only: replace the largest upper bound with +inf and the smallest negative lower bound with -inf before sampling, so a reaction whose biological maximum exceeds the model's arbitrary cap is not pinned there. Off by default.

min_flux:

random_objective only: re-solve each sample parsimoniously to minimise total flux at the optimum, which removes residual loop flux.

n_proc:

random_objective only: worker processes for drawing samples in parallel (each sample is an independent LP solve). Defaults to cobra.Configuration().processes; 1 draws the samples serially in this process. The result is reproducible for a given seed regardless of n_proc: each sample draws from its own independent RNG stream spawned from seed.

Returns

FluxSamplingResult :

reporter_metabolites

reporter_metabolites(model: cobra.Model, gene_pvalues: Mapping[str, float], *, gene_fold_changes: Mapping[str, float] | None=None)

Compute Reporter Metabolites from per-gene differential-expression p-values.

gene_pvalues maps gene id to p-value. Genes not in the model and genes with a NaN (or None) p-value are dropped. Any other p-value outside [0, 1], including ±inf and for genes that are not in the model, raises ValueError (RAVEN raises RAVEN:badInput). If gene_fold_changes (gene id to log fold change) has at least one non-zero value, two extra results are returned for the up- (fc ≥ 0) and down- (fc < 0) regulated gene subsets, in addition to "all"; with all fold changes zero (or none given) only "all" is returned, as in RAVEN. A metabolite with a single scored neighbouring gene has std_z 0.

The background correction uses the population standard deviation (ddof=0) of the scored Z-scores, which the closed form requires. The reported p_value is the one-sided enrichment 1 - Φ(z) and each table is sorted by z_score in descending order.

Parity with RAVEN's reporterMetabolites: the same closed-form correction, one-sided p-value and sort order, with "up" and "down" matching its "only up" and "only down" entries.

Parameters

model : cobra.Model

Model whose metabolites are scored through the genes of their neighbouring reactions.

gene_pvalues : Mapping[str, float]

Gene id to differential-expression p-value, as described above.

gene_fold_changes : Mapping[str, float] or None

Gene id to log fold change. The default None returns only the "all" result. A gene with a p-value but no entry here, or a NaN entry, is in neither subset (RAVEN takes one fold change per gene, and a NaN fold change satisfies neither >= 0 nor < 0).

ReporterResult

Reporter-metabolite scores for one gene set.

test is "all", "up" or "down"; table is a DataFrame with columns metabolite, name, z_score, p_value, n_genes, mean_z, std_z sorted by descending z_score.

trace_flux_path

trace_flux_path(model: cobra.Model, fluxes: pd.Series | Mapping[str, float], from_rxn: str, to_rxn: str, *, cutoff: float=_DEFAULT_CUTOFF, max_hops: int=_DEFAULT_MAX_HOPS, trace_material: bool=True, carbon_only: bool=False, exclude_mets: Iterable[str] | None=None)

Find the highest-flux-fraction path from from_rxn to to_rxn.

Parameters

model:

Model both reactions belong to.

fluxes:

Flux for every reaction the search might visit, indexed by reaction id (a full-model flux vector, for example from model.optimize().fluxes, always qualifies).

from_rxn, to_rxn:

Source and target reaction ids.

cutoff:

A reaction counts as consuming a metabolite when its stoichiometric coefficient times its flux is below -cutoff.

max_hops:

Longer paths are pruned.

trace_material:

Exclude common currency and cofactor metabolites (ATP, ADP, NAD, NADH, H2O, CoA, Pi, CO2, ...) so the path follows material flow.

carbon_only:

Additionally require each intermediate metabolite to contain at least one carbon atom (via Metabolite.elements, which parses the formula, so for example a calcium ion is not counted as carbon-bearing). A metabolite with no formula counts as having no carbon.

exclude_mets:

Extra metabolite names or ids (compartment suffix stripped) to exclude, on top of the currency list.

Returns

TraceFluxPathResult :

Raises

ValueError

If from_rxn or to_rxn is not in the model, or fluxes has no entry for a reaction the search needs.

TraceFluxPathResult

Outcome of a trace_flux_path search.

Parameters

reactions:

Reaction ids on the best path, including from_rxn and to_rxn. Empty if no path was found.

metabolites:

Metabolite id at each junction; len(metabolites) == len(reactions) - 1.

cumulative_fraction:

Fraction of from_rxn's flux that passes through this exact path (0 = no path found, 1 = all of it reaches to_rxn directly).

walk_fluxes

walk_fluxes(model: cobra.Model, fluxes: pd.Series | Mapping[str, float], start_rxn: str | int, *, cutoff: float=1e-08, max_per_met: int=8, input_fn=input, print_fn=print)

Interactively navigate a flux distribution reaction by reaction.

Starting from start_rxn, prints all flux-carrying reactions connected through shared metabolites, grouped by metabolite and labelled with each neighbour's role (produces or consumes). Entering a neighbour's number steps to it, b goes back, and q quits. input_fn and print_fn are injectable for scripting or testing; they default to the builtins.

See FluxWalker for a non-blocking, programmatic form of the same navigation (.step(n), .back(), .groups).

Parameters

model : cobra.Model

Model whose reactions are navigated.

fluxes : pd.Series or Mapping[str, float]

Flux by reaction id. A reaction without an entry counts as carrying zero flux.

start_rxn : str or int

Reaction to start from, as a reaction id or as an index into model.reactions. An id that is not in the model raises ValueError.

cutoff : float

A metabolite is shown when the absolute net contribution of the current reaction is at least this value and at least one neighbour has an absolute contribution above it. A neighbour is listed when its absolute contribution exceeds it. The default is 1e-8.

max_per_met : int

Maximum number of neighbouring reactions listed per metabolite, taken in order of descending absolute flux. The default is 8.

input_fn : callable

Function called with the prompt "  > " that returns the user's response as a string. The default is the builtin input.

print_fn : callable

Function called with each output string, which can contain line breaks. The default is the builtin print.