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: carry real (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 — entry point for all 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, choosing the sign by hand (the p-value alone 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 flux and expression both changed significantly, in the same direction (RAVEN's pR).

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

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

Notes

Two simplifications relative to analyzeSampling.m: it computes an "opposite flux direction" branch (pM=erf(|Z|), pH=(1-pM)*..., pR=0) whose results are unconditionally overwritten by the very next branch before scores is ever read — dead code, verified by tracing every assignment, not replicated here. And it bundles its own tcdf (a pure-MATLAB fallback for when the Statistics Toolbox isn't installed); this uses scipy.stats.t.cdf directly, which is exact and always available.

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 (e.g. wild type), indexed by reaction id. Every reaction in the model must be present.

comparison:

Flux for the condition being compared against it (e.g. 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 left out of the result entirely rather than reported as unchanged. A name matching no metabolite is reported in CompareFluxesResult.missing_metabolites rather than raising.

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 the reference flux is below the cutoff — the reaction has no baseline to be relative to. type is "on", "off", "flip" or "".

turned_on / turned_off / flipped are reaction-id lists in the model's own order, not the sorted order of changed.

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: carry real (non-loop) flux.

A reaction is kept if its FVA range spans more than 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 — the right test for "loopy", unlike a fixed bound threshold which wrongly drops legitimate reactions that simply reach the model's default (e.g. 1000) bound. Pass loopless=False for a faster, looser pass that keeps any flux-carrying reaction (loops included).

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. A False here is not fatal (the last ellipsoid iterate is still a valid rounding) but very elongated results may warrant more thinning.

n_warmup:

achr only — number of FVA warmup directions.

fixed_reactions:

chrr only — reactions folded into the equality system as fixed.

FluxWalker

Steppable navigator over a solved flux distribution.

fluxes is indexed by reaction id (e.g. 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)

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

Enforces target flux from max_fraction/n_steps up to max_fraction of the theoretical maximum in n_steps steps, maximising growth (biomass_rxn or the model's current objective) with pFBA at each step. Returns an FSEOFResult.

Parameters

correlation_threshold:

Minimum absolute correlation between a reaction's flux and the enforced product flux for it to be reported as a target; raise it to keep only the most cleanly-trending reactions.

flux_eps:

Tolerance below which a flux or slope is treated as zero.

min_target:

Floor of the enforced-flux range, evenly spaced up to the scan's ceiling (max_fraction of the theoretical maximum) in n_steps steps. Defaults to None, which floors the range at max_fraction of the theoretical maximum / n_steps instead. Pass e.g. the target reaction's flux at growth-optimum to scan from there.

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 (e.g. 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 gets exactly ±100 rather than a division by zero.

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; constrain it yourself first — e.g. pin it near its optimum — if a minimal solution among optimal ones is wanted).

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 (more negative = more discouraged); non-negative entries are clamped up to the largest (least negative) negative score present, matching MATLAB's "positive scores are not possible" rule — or to 0 if scores has no negative entries at all, as RAVEN's getMinNrFluxes does. Defaults to an implicit -1 for every reaction (plain unweighted cardinality minimization).

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 that keeps the two Newton solves well-conditioned.

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/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 (e.g. "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.

A simpler, Pythonic rendering of traceFluxPath.m's console diagram — not a literal reproduction of its text; the structured result is what's compared for parity.

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 — entry point for all 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/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"/"chrr" methods draw the (near-)uniform interior distribution; "random_objective" draws diverse vertices. Set any constraints you want to condition on (e.g. a biomass lower bound, measured fluxes, enzyme-usage bounds) on the model before calling.

Parameters

n_samples:

Number of flux vectors to return.

method:

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

seed:

Seed for reproducible chains/draws.

thinning:

achr/chrr — Markov-chain steps between recorded samples (default 100). The default is calibrated for small models (~200–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 roughly 12 effective samples from 300 stored samples (ESS ≈ n × (1−ρ)/(1+ρ)). Independent chains at these default settings also frequently fail to agree with each other: on yeast-GEM, the median reaction fails the Gelman-Rubin R-hat>1.1 convergence threshold across 4 independent chains (67.5% of reactions fail it; 96.5% fail the stricter R-hat>1.01). Treat default-settings genome-scale output as unconverged for most reactions, not a caveat for a minority of hard ones. For genome-scale analyses either increase thinning substantially (≥1000), increase n_samples to compensate, try method='chrr' (raven-toolbox's rounding-based alternative, not wired into ACHR), or use an ESS/R-hat diagnostic to assess sample quality (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 folded into the equality system as fixed (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 finite bound with infinity before sampling, so a reaction whose biological maximum exceeds the model's arbitrary cap isn't pinned there. Off by default.

min_flux:

random_objective only — re-solve each sample parsimoniously to minimise total flux at the optimum, squeezing residual loop flux out.

n_proc:

random_objective only — worker processes for drawing samples in parallel (each sample is an independent LP solve). Defaults to cobra.Configuration().processes; set to 1 to sample serially. Reproducible for a given seed regardless of n_proc: each sample draws from its own independent RNG stream spawned from seed, not from one continuously-advancing generator.

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 → p-value (genes not in the model, or with a NaN or out-of-[0, 1] p-value, are dropped — a stray invalid p-value would otherwise turn the whole result NaN). If gene_fold_changes (gene id → log fold change) is given, two extra results are returned for the up- (fc ≥ 0) and down- (fc < 0) regulated gene subsets, in addition to "all".

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

Parity with RAVEN's reporterMetabolites: the same closed-form correction, one-sided p-value and sort order, with "up"/"down" matching its "only up"/"only down" entries. A RAVEN version that divides by the sample standard deviation (std with N−1) gives z-scores smaller by a factor √(N/(N−1)) for N scored genes; the ranking is the same.

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, e.g. from model.optimize().fluxes, always qualifies).

from_rxn, to_rxn:

Source and target reaction ids.

cutoff:

Minimum |flux| for a reaction to count as consuming a metabolite.

max_hops:

Longer paths are pruned.

trace_material:

Exclude common currency/cofactor metabolites (ATP, ADP, NAD, NADH, H2O, CoA, Pi, CO2, ...) so the path follows material flow rather than an energy-carrier shortcut.

carbon_only:

Additionally require each intermediate metabolite to contain at least one carbon atom (via Metabolite.elements, which parses the formula properly — unlike a naive regex, it does not mistake e.g. a calcium ion for a carbon-bearing one). A metabolite with no formula counts as having none.

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 routes 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/consumes). Enter a neighbour's number to step to it, b to go back, or q to quit. input_fn/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).