diff --git a/src/raven_toolbox/analysis/sampling.py b/src/raven_toolbox/analysis/sampling.py index cbb123c..c3dbaa7 100644 --- a/src/raven_toolbox/analysis/sampling.py +++ b/src/raven_toolbox/analysis/sampling.py @@ -104,18 +104,10 @@ def random_sampling( 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 + On genome-scale models (>2000 reactions) consecutive samples remain strongly + correlated at the default, and independent chains can disagree. Options for + genome-scale analyses are increasing thinning substantially (≥1000), increasing + n_samples, or assessing sample quality with an ESS or R-hat diagnostic (see the `sampling convergence calibration study `_ on raven-docs). diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 8c5ebfc..6bcc8c0 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -97,12 +97,8 @@ def _reaction_big_m(lb: float, ub: float, cap: float, force_on: float) -> float: # A resolve_ties tie-break phase is adopted whether it proves or times out (see # _resolve_ties), so it is given a fraction of the primary solve's own budget rather than -# the full amount again: a phase that would not converge within the full budget has been -# observed (on genome-scale Human-GEM) to still not converge, only slower; a phase that -# does converge does so far inside either budget (a pure integer count/id-rank objective -# with a 0.4 absolute gap). Not independently tuned — see the ftINIT reproducibility study -# on raven-docs for the experiment that would calibrate this properly instead of halving -# on general principle. +# the full amount again: a phase that does converge does so far inside either budget (a +# pure integer count/id-rank objective with a 0.4 absolute gap). Not tuned. _TIE_BREAK_TIME_FRACTION = 0.5 @@ -223,15 +219,12 @@ def run_ftinit( prove_abs_gap : float or None When set, proves the solve to this fixed absolute objective gap (on Gurobi: ``MIPGap`` 0 and ``MIPGapAbs`` this value); ``mip_gap`` and ``mip_gap_abs`` do not apply then. The - default None leaves the gap to ``mip_gap`` and ``mip_gap_abs``. In the four-seed sweep - (seeds 1234, 7, 1 and 2; ``resolve_ties=True``; no gap-fill; Human-GEM v2.0.0-38, one - cell line) the default escalation returned a step-1 objective 4.0 to 12.0 below the - proven optimum without a warning, while every value from 2.0 to 0.05 reached the - optimum for all four seeds and gave the same reaction sets. Every step solve (the 48 - primary solves) ended with status ``optimal``; the flux phases ``phase2c-minflux`` (3 - of 4 seeds) and ``phase2d-direction`` (4 of 4 seeds) reached the time limit and were - adopted as unproven incumbents. The docstring of :func:`raven_toolbox.init.ftinit` - gives the measured run times. + default None leaves the gap to ``mip_gap`` and ``mip_gap_abs``. The default escalation + can stop at an objective below the proven optimum without a warning. The flux phases + ``phase2c-minflux`` and ``phase2d-direction`` can reach the time limit, in which case + the incumbent is adopted as an unproven tie-break. The docstring of + :func:`raven_toolbox.init.ftinit` links the study that measures the effect of this + setting. resolve_ties : bool Opt-in; the default False gives RAVEN's result. When True, the degeneracy of the MILP is resolved deterministically. After the score optimum is found, a lexicographic @@ -282,11 +275,10 @@ def run_ftinit( Gurobi's ``Threads`` (default 1, RAVEN's value). More threads can shorten the solve, and with more than one thread two runs of the same input can return different reaction sets, because multi-threaded Gurobi chooses between equal optima non-deterministically. - The studies on raven-docs ran with one thread. The number of threads has no measured - effect on the spread between different seeds, so it does not replace ``resolve_ties``; - see the `ftINIT reproducibility study + The number of threads does not select between equal optima in a controlled way, so it + does not replace ``resolve_ties``. The `ftINIT reproducibility study `_ - on raven-docs. It is applied on Gurobi only. + on raven-docs measures the spread between seeds. It is applied on Gurobi only. feas_tol : float Gurobi's ``FeasibilityTol``. The default is 1e-9 (RAVEN's value and Gurobi's minimum; Gurobi's own default is 1e-6). See the note on tolerances below. @@ -305,8 +297,7 @@ def run_ftinit( ``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. + run time. The private parameter ``_pin_fluxes`` (default True) is set by :func:`raven_toolbox.init.ftinit`: it is False in the last step, whose fluxes no later @@ -635,9 +626,7 @@ def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit, rank_key Many reaction subsets reach the same score optimum of the ftINIT MILP, and the solver returns one of them. The choice is reproducible for a fixed solver build and seed, and can - change with the Gurobi version, the platform or the seed. In the determinism study - (Human-GEM v2.0.0-38, one cell line) two seeds gave kept-reaction sets that differ by 110 - reactions and 17 predicted essential genes. + change with the Gurobi version, the platform or the seed. The function runs a lexicographic phase 2 that holds the score objective at its optimum with a floor constraint: first it minimises the count of kept removable reactions (the @@ -655,14 +644,11 @@ def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit, rank_key distribution, which the staged schedule reads for the forced directions and the ``ess_force`` clamp of the next step, is pinned afterwards by :func:`_min_abs_flux`. - Caveat, measured on genome-scale Human-GEM (both phases prove in milliseconds at toy or - unit-test scale): the phases can exhaust ``time_limit``, in which case the incumbent - adopted is an unproven tie-break. It reduces the run-to-run spread relative to no - tie-break at all, and it does not guarantee a unique result. A warning is issued when this - happens; see the `ftINIT reproducibility - study + Caveat: the phases can exhaust ``time_limit``, in which case the incumbent adopted is an + unproven tie-break. It does not guarantee a unique result. A warning is issued when this + happens. The `ftINIT reproducibility study `_ - on raven-docs. + on raven-docs measures the spread between seeds with and without the tie-break. Returns ``True`` if phase 2 produced a usable solution (the caller then reads the resolved "on" set and fluxes), ``False`` if it did not converge to an incumbent (the @@ -1158,21 +1144,12 @@ def ftinit( relative-gap escalation with a single solve per step proven to this fixed absolute gap, so the kept-set objective is optimal at any scale and not an arbitrary within-gap incumbent. It changes which optimum is returned and not its score; the alternative - optima are equally consistent with the expression data. The sweep ran on Human-GEM - v2.0.0-38 with one cell line, four seeds (1234, 7, 1 and 2), ``resolve_ties=True``, - no gap-fill and ``time_limit`` 1800 s. The default escalation returned a step-1 - objective 4.0 to 12.0 below the proven optimum; with seeds 1234 and 7 alone it ended - 7.0 and 12.0 below, and 4.0 is seed 2. Every value from 2.0 to 0.05 reached the - optimum for all four seeds and gave the same reaction sets, so 1.0 lies in the range - of equal results; 5.0 stopped 2.0 below the optimum for one seed. Every step solve - (the 48 primary solves) ended with status ``optimal``; the flux phases - ``phase2c-minflux`` (3 of 4 seeds) and ``phase2d-direction`` (4 of 4 seeds) reached - the time limit and were adopted as unproven incumbents. The mean build time was about - 5770 s at 0.05 against 5640 s at 1.0. The study recommends 1.0 together with - ``resolve_ties`` when the seed-to-seed swing of ``resolve_ties`` alone (41 reactions, - 1 gene here) is too large and the run time is acceptable. The full-pipeline runs of - the calibration study used ``time_limit`` 600 s per step and the determinism builds - 1800 s. + optima are equally consistent with the expression data. A larger value stops the + solve at a larger distance from the proven optimum. The flux phases of ``resolve_ties`` + can reach ``time_limit`` and are then adopted as unproven incumbents. The `ftINIT + reproducibility study + `_ + on raven-docs sweeps this setting and reports the run times. resolve_ties : bool Opt-in; the default False gives RAVEN's behaviour. It reduces the arbitrariness of the tie-break in the degenerate MILP and gives a more parsimonious and more reproducible @@ -1206,20 +1183,13 @@ def ftinit( Notes ----- - The studies on raven-docs ran at the default tolerances, ``presolve`` and one thread (see - :func:`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 + The stability of the extraction under small input changes is a separate property that + neither ``seed`` nor ``threads`` determines. The `ftINIT reproducibility study `_ - 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. + on raven-docs measures the spread between seeds and the run time for the default + escalation, for ``resolve_ties`` and for ``resolve_ties`` with ``prove_abs_gap``. The + tie-break phases of ``resolve_ties`` can exhaust ``time_limit``; the result is then an + unproven tie-break, not a unique optimum, 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 diff --git a/src/raven_toolbox/init/prep.py b/src/raven_toolbox/init/prep.py index b943d5a..ed2e17f 100644 --- a/src/raven_toolbox/init/prep.py +++ b/src/raven_toolbox/init/prep.py @@ -323,9 +323,10 @@ def prep_init_model( larger merged model that may reach a different extraction optimum. scale : bool Keyword-only; default ``True``. When ``False``, the merged model's stoichiometry is - not rescaled. In the INIT parameter calibration study on raven-docs (Human-GEM, one - cell line), the full ftINIT pipeline without scaling was feasible at ``big_m`` 1000 - and took 157 s against 83 s with scaling, and it was infeasible at ``big_m`` 100. + not rescaled. Without scaling, a small ``big_m`` can make the MILP infeasible. The + `INIT parameter calibration study + `_ + on raven-docs compares the pipeline with and without scaling. verbose : bool Keyword-only; default ``False``. Prints one line per phase plus per-task progress during essential-reaction discovery. diff --git a/src/raven_toolbox/init/taskfill.py b/src/raven_toolbox/init/taskfill.py index 852df0e..3127ad4 100644 --- a/src/raven_toolbox/init/taskfill.py +++ b/src/raven_toolbox/init/taskfill.py @@ -41,12 +41,8 @@ # A resolve_ties tie-break phase is adopted whether it proves or times out (see # _resolve_ties_fill), so it is given a fraction of the primary fill's own budget rather -# than the full amount again: a phase that would not converge within the full budget has -# been observed (on genome-scale Human-GEM) to still not converge, only slower; a phase -# that does converge does so far inside either budget (a pure integer count/id-rank -# objective with a 0.4 absolute gap). Not independently tuned — see the ftINIT -# reproducibility study on raven-docs for the experiment that would calibrate this -# properly instead of halving on general principle. +# than the full amount again: a phase that does converge does so far inside either budget +# (a pure integer count/id-rank objective with a 0.4 absolute gap). Not tuned. _TIE_BREAK_TIME_FRACTION = 0.5 diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index 4f45a8e..724886e 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -131,13 +131,11 @@ def _prepare_scope(model, scores, reactions_to_relocate, *, transportable, base_ # --------------------------------------------------------------------------- master def _pin_deterministic(prob, opt) -> None: """Pin solver parameters so a degenerate placement objective gives the same solution on every - run and on every machine. In a benchmark of the yeast master (identical model, one parameter - varied at a time), three settings changed which co-optimal placement Gurobi returned: thread - count, seed and presolve level. The MIP gap changed the result only when loosened above ~1e-2, - which also gives non-optimal solutions, and no tolerance changed it. The function sets one - thread, a fixed seed, presolve level 2 (the optimizeProb default of the MATLAB port), a MIP gap - of 0 (exact optimum) and an integrality tolerance of 1e-9. Only the Gurobi backend is - configured; other backends keep their defaults.""" + run and on every machine. The objective is degenerate, so thread count, seed and presolve level + can change which co-optimal placement Gurobi returns. The function sets one thread, a fixed + seed, presolve level 2 (the optimizeProb default of the MATLAB port), a MIP gap of 0 (exact + optimum) and an integrality tolerance of 1e-9. Only the Gurobi backend is configured; other + backends keep their defaults.""" if "gurobi" not in getattr(prob, "__name__", ""): return try: diff --git a/src/raven_toolbox/localization/scores.py b/src/raven_toolbox/localization/scores.py index add37f1..c57a285 100644 --- a/src/raven_toolbox/localization/scores.py +++ b/src/raven_toolbox/localization/scores.py @@ -193,8 +193,7 @@ def load_deeploc(path: str | Path, *, labels). A label absent from the map is dropped with its column (see the module docstring). min_confidence : float Drops genes whose top compartment probability (before per-gene normalisation) is below it - (default ``0.0``, which keeps every gene). In the yeast-GEM benchmark of the confidence - study, the share of calls that other evidence corroborates rises with the raw probability. + (default ``0.0``, which keeps every gene). membrane_split : mapping of str to str or None For example ``{"m": "mm"}`` (default ``None``, no split). It moves the probability of an organelle to its membrane sub-compartment when the protein is membrane-associated @@ -202,7 +201,7 @@ def load_deeploc(path: str | Path, *, it stays in the lumen. Keys and values are the compartment ids of the model (after ``compartment_map``). A key that is not a compartment of the table is skipped. Raises :class:`ValueError` when the file has no ``Soluble`` column. The evidence supports only the - mitochondrial split (``m``/``mm``; matrix versus membrane AUC ~0.93). DeepLoc does not + mitochondrial split (``m``/``mm``). DeepLoc does not separate the ER lumen from the membrane, so a split of ``er`` has no support. membrane_threshold : float Value of ``1 - P(Soluble)`` from which a protein counts as membrane-associated (default @@ -214,9 +213,8 @@ def load_deeploc(path: str | Path, *, normalise : bool If True (default), rescale the best compartment of every gene to 1.0. This follows the RAVEN ``parseScores`` convention and keeps the scales of several sources comparable. False - keeps DeepLoc's raw probabilities. The raw probabilities are calibrated (a 0.97 call is - more reliable than a 0.40 one), so a downstream assignment can weight confident genes more - heavily. + keeps DeepLoc's raw probabilities, so a downstream assignment can weight confident genes + more heavily. """ raw = pd.read_csv(path) wide = _wide_mapped(raw, id_column=None, compartment_map=compartment_map, source=str(path)) diff --git a/src/raven_toolbox/localization/triage.py b/src/raven_toolbox/localization/triage.py index cf7c679..68e9b7f 100644 --- a/src/raven_toolbox/localization/triage.py +++ b/src/raven_toolbox/localization/triage.py @@ -8,9 +8,7 @@ ``uncertainty`` in ``[0, 1]`` blended from the signals that apply, and the list is ranked by it. The signals that fire differ per gene: -* **low confidence**: the raw top probability of the predictor is low (in the yeast-GEM - benchmark, the share of calls that other evidence corroborates rises with the raw - probability). This signal carries the largest default weight (0.35) of the signals in the +* **low confidence**: the raw top probability of the predictor is low. This signal carries the largest default weight (0.35) of the signals in the blended ``uncertainty``. Needs ``scores.raw_confidence`` (``load_deeploc(keep_raw_confidence=True)``) or an explicit ``raw_confidence=``, because normalisation hides the raw probability. @@ -28,7 +26,8 @@ See the `DeepLoc yeast-GEM benchmark `_ -on raven-gecko-parity for the source of the thresholds. +on raven-gecko-parity compares the raw predictor probability with other evidence for the +yeast-GEM calls and gives the bins and thresholds. """ from __future__ import annotations @@ -61,7 +60,7 @@ def confidence_bin(p: float) -> str: - """Calibration bin of a raw confidence, with the bins of the benchmark's corroboration rates.""" + """Bin of a raw confidence, with the bins of the benchmark's corroboration table.""" for lo, hi in ((0.0, 0.5), (0.5, 0.7), (0.7, 0.9)): if lo <= p < hi: return f"[{lo:.1f},{hi:.1f})" diff --git a/src/raven_toolbox/reconstruction/homology/homology.py b/src/raven_toolbox/reconstruction/homology/homology.py index 2108d40..2239f16 100644 --- a/src/raven_toolbox/reconstruction/homology/homology.py +++ b/src/raven_toolbox/reconstruction/homology/homology.py @@ -352,14 +352,10 @@ def get_model_from_homology( exclude every hit that it would exclude. min_align_len : int A hit passes only when its alignment length is at least this value. Default 100, as - RAVEN's ``minLen``. Values at or below 150 measured the same; on each of the four - organisms tested, the score against KEGG annotations (wrong matches counted double) - was 0.010 to 0.020 lower at 200 than at 100. + RAVEN's ``minLen``. min_identity : float - A hit passes only when its identity is at least this value. Default 40. The value - lies in the 35 to 45 range that the two reference sources favour; 35 scores higher - for two organisms (``yli`` 0.868 against 0.860, ``ani`` 0.814 against 0.804). It is - the filter that limits the result. + A hit passes only when its identity is at least this value. Default 40. It is the + filter that limits the result. review_identity : float or None Lower identity threshold (for example 25) that collects the reactions the thresholds reject. Any reaction that transfers at this lower identity but not at diff --git a/src/raven_toolbox/reconstruction/kegg/query.py b/src/raven_toolbox/reconstruction/kegg/query.py index 05859de..69cb1c4 100644 --- a/src/raven_toolbox/reconstruction/kegg/query.py +++ b/src/raven_toolbox/reconstruction/kegg/query.py @@ -56,8 +56,8 @@ def run_hmmsearch( One ``hmmsearch`` of the whole concatenated multi-profile ``library`` (query) against ``fasta`` (target), which is the fast search direction. It runs with ``--cpu`` and needs no ``hmmpress``. ``-Z`` is set to the profile count (when it - is positive) so that the per-hit E-values match the convention :func:`assign_kos` is - calibrated against (identical to a ``hmmscan`` against the same library). A + is positive) so that the per-hit E-values match the convention that :func:`assign_kos` + expects (identical to a ``hmmscan`` against the same library). A ``RuntimeError`` is raised if ``hmmsearch`` exits with a non-zero status. Parameters @@ -134,21 +134,13 @@ def assign_kos( negative, the best (smallest) hit has ratio 1, and weaker hits have a smaller positive ratio. - The defaults were calibrated by cross-validation against the KEGG gene and KO - annotation of four organisms that span the prokaryote and eukaryote libraries and - well- and less-studied species (*S. cerevisiae*, *Cyanidioschyzon merolae*, - *E. coli*, *Mycoplasma genitalium*). Annotated genes have a median E-value of - 1e-100 to 1e-155, and spurious hits cluster at about 1e-8, a separation of ~20 - orders of magnitude. RAVEN's cut-off of ``1e-50`` lies within the E-value range of - the annotated genes and drops real but divergent hits. For the divergent minimal genome - (*M. genitalium*, ``mge``) the gene-to-KO recall is 0.95 at 1e-30 and 0.84 at 1e-50, - with no gain in rejecting noise (noise E-values are far weaker). The default cut-off - is therefore 1e-30, which recovers those hits and is ~22 orders of magnitude above - the noise level. ``min_score_ratio_g`` (default 0.9) sets the precision, because it - resolves genes with several KOs. Against RAVEN's settings (1e-50 and 0.8), the - recall for ``mge`` rises from 0.84 to 0.94. ``min_score_ratio_ko`` values of 0.0, 0.3 - and 0.5 changed precision and recall by at most 0.02 across the four organisms; 0.3 is - kept as in RAVEN. + The default cut-off is 1e-30, weaker than RAVEN's ``1e-50``, so that real but divergent + hits are kept; the ratio filters then remove weak members. ``min_score_ratio_g`` (default + 0.9, RAVEN uses 0.8) sets the precision, because it resolves genes with several KOs. + ``min_score_ratio_ko`` (default 0.3) is kept as in RAVEN. + The `KEGG HMM cut-off calibration study + `_ + on raven-docs compares these settings with the KEGG gene and KO annotation of four organisms. Parameters ----------