feat: add human_gem Python package with load_model() - #1069
Conversation
First step towards a pip-installable Python interface to Human-GEM, on
raven-toolbox + cobrapy. human_gem.load_model() loads Human-GEM as a fully
annotated cobra model - reading the YAML via raven-toolbox and merging the
cross-references and SBO terms from the annotation tables, which a bare
cobra.io.load_yaml_model does not do:
import human_gem
model = human_gem.load_model() # 12877 rxns, 8460 mets, cross-refs attached
- pyproject.toml (deps: cobra, pandas, raven-toolbox)
- human_gem/io.py load_model(model_dir=None, annotate=True)
- human_gem/annotation.py TSV -> MIRIAM/SBO merge (packaged home of the logic
currently duplicated in code/annotateGEM.py)
Verified locally: load_model() attaches cross-references to 12698/12877
reactions and 8449/8460 metabolites; annotate=False returns the bare model.
Next increments: reduce code/annotateGEM.py to a shim importing from human_gem
(single source of truth) + install the package in CI; then promote the ftINIT
essentiality and task helpers into the package.
Model quality report
Each check name links to its explanation in the testResults README. Model checksDuplicate keys (model unloadable) and no growth block the merge; every other row is a non-blocking report.
MACAW and mass/charge balance
Model file and metabolic tasks
MEMOTETotal score: 63.2% (core subset) 0
Per-test scores
Full suite: 64.2% 0 · from the last The score above is the fast core subset. Comment Gene essentiality (Hart 2015)Not run automatically (it takes hours). Comment ❌ = a count rose vs the target branch (regression) · Full workflow run · this comment is edited as results come in |
raven-toolbox 3.0.0b1 is on PyPI, so the package no longer needs a git URL. The specifier names the pre-release because pip will not resolve one otherwise. The same release adds two ftINIT settings that the gene-essentiality run now uses. The extraction MILP is degenerate, so which of the equally-scoring reaction sets comes back was left to the solver, and the predicted essential genes moved with it: a develop-versus-branch comparison showed most changed genes flipping in anticorrelated directions across cell lines, unrelated to the edit under review. resolve_ties takes the sparsest, then lowest-id, optimum, and prove_abs_gap=1.0 proves each step to a fixed absolute gap rather than accepting an incumbent within a relative one. Neither makes the model biologically more accurate; both make a re-run report the same genes.
The two regenerated result files conflicted; taken from develop so the graded gene-essentiality documentation added by #1077 is kept.
|
Gene essentiality was run on this pull request via GH Actions. Essentiality for a viability task ( Graded gene essentiality vs Hart 2015 (task-scope analysis)
The full per-gene matrix is committed to
|
…1079) * feat: parallelize gene-essentiality across cell lines via a build matrix The gene-essentiality workflow built and scored all five Hart 2015 cell-line models in one sequential loop in a single job, taking hours end to end even though every cell line is independent of the others. `gradedEssentiality.py`: had no way to build a subset of cell lines, so splitting the work across parallel processes or jobs was not possible. Adds `--shard-index`/`--shard-count` to select a slice of cell lines, and `--aggregate-only` to combine every cell line's checkpoint into the final report without rebuilding any of them. `gene-essentiality.yml`: ran Step 1 (prepINITModel) and every cell line's ftINIT build and essentiality scan in one job, so wall-clock time scaled with the number of cell lines. Splits into a `build` job matrixed one shard per cell line and an `aggregate` job that combines their checkpoints into the committed report; the `diagnose` path is unchanged but now runs as its own job instead of an early exit inside the shared script. * chore: update model QC results [skip ci] * perf: dedupe and parallelize the gene-essentiality knockout scan The per-tissue scan in taskEssentialGenes.py solved thousands of small, nearly identical LPs sequentially, so its runtime did not benefit from genes sharing a knockout or from the extra cores a CI runner has. - `_scan`: genes that disable the identical reaction set (subunits of the same AND-linked complex, e.g. an ETC complex) were each solved separately against every task instead of once per distinct reaction set. - `_task_feasible_without`'s solves ran under Gurobi's default automatic/ concurrent method, which spins up several threads per solve for a problem this small, adding synchronisation overhead across thousands of solves for no benefit. - `find_task_essential_genes`/`find_task_essential_categories` had no way to spread the scan's independent gene knockouts across more than one core. Verified against the previous per-gene scan (sequential and processes=2) on a synthetic model covering a lone-gene reaction, an isoenzyme pair and a multi-subunit complex. * fix: run the knockout-scan pool with a spawn context `_scan`'s `ProcessPoolExecutor` used the platform default multiprocessing context, which is `fork` on the Linux CI runners this targets. Forking after `base`'s Gurobi environment is already running duplicates that environment, and any internal license/logging threads it holds, into every worker; this is the same non-fork-safe state that made cobra's default parallel FVA deadlock on Linux CI runners, which is why `estimateEssentialGenes.estimate_essential_genes` and `gradedEssentiality.main` both pin FVA to `processes=1`. - `_scan`: pass `mp_context=multiprocessing.get_context("spawn")` to the pool so each worker starts as a fresh interpreter and only creates its own Gurobi environment inside `_init_worker`, after the fork boundary, leaving nothing shared between workers to deadlock on. * perf: scan each cell line's gene knockouts on all of the runner's cores `gradedEssentiality.py` shards cell lines across CI jobs (#1079), but within one cell line's `find_task_essential_categories` call every gene knockout still ran on a single core, leaving a shard's other cores idle regardless of how many the runner had. - `main`: adds `--processes` (default: `os.cpu_count()`) and forwards it to `find_task_essential_categories`, which already supports splitting its independent gene knockouts across a process pool. * chore: update model QC results [skip ci]
The gene-essentiality run at https://github.com/SysBioChalmers/Human-GEM/actions/runs/34445477424 computed results for all five cell lines (posted as a PR comment) but its "Auto-commit results" step failed to push: "Updates were rejected because the tip of your current branch is behind its remote counterpart". - `gene-essentiality.yml`: "Update local branch before committing changes" ran `git fetch` then `git checkout $BRANCH_NAME`, but `git fetch` only moves the `origin/$BRANCH_NAME` ref; the local branch stayed at whatever commit the job started from. Anything pushed to the branch while the job ran left the auto-commit step basing its commit on that stale tip, so its push was rejected as non-fast-forward and the finished run's results were lost.
* chore: restore generated model files removed by #1061 * chore: new version * chore: add the release workflow so it can be dispatched workflow_dispatch resolves against the default branch, so a Release workflow that exists only on develop cannot be started. Adding it here makes it available in the Actions tab. The workflow calls code/io/increaseHumanGEMVersion.py from the develop checkout it cuts the release branch from, so a release can only run once #1078 is on develop. * ci: edit the /run gene-essentiality comment instead of replying `/run gene-essentiality` posted its "started" acknowledgement as a separate new comment, so a PR thread that used the command showed two comments (the command and the acknowledgement) instead of one. - `pr-command.yml`: the `gene-essentiality` job's "started" notice is now appended to the triggering comment itself (`gh api ... -X PATCH`) rather than posted with `gh pr comment`.
Each of the 5 shards in `main` independently called `_prep_human_model_for_ftinit` whenever it had a pending cell line, even though that step's output does not depend on which cell line asks for it (it runs before any cell line's expression data is applied), so the same tens-of-minutes computation ran five times over. - `gradedEssentiality.py`: adds `_load_or_build_prep`, `--prep-cache` and `--prep-only`. `--prep-cache PATH` loads Step 1's `PrepData` from `PATH` if present instead of rebuilding it, and writes it there otherwise; `--prep-only` builds it and exits without touching any cell line. Without `--prep-cache`, a run still builds its own Step 1, as before.
The `build` matrix ran `gradedEssentiality.py` for each of its 5 shards without sharing anything between them, so the newly cacheable Step 1 still ran once per shard by default. - `gene-essentiality.yml`: adds a `prep` job that builds `PrepData` once with `--prep-only --prep-cache` and uploads it as an artifact; `build` now depends on `prep`, downloads that artifact, and passes `--prep-cache` so every shard loads the shared prepData instead of rebuilding its own.
|
/run gene-essentiality |
|
🧬 Started gene-essentiality on |
The full metabolic task list was missing the leading (empty) column that metabolicTasks_Essential.txt has, so RAVEN parseTaskList read the whole file as a single task. Add the leading column to align the format; parseTaskList now returns 257 tasks and compareMultipleModels runs against the full list. No task definitions are changed.
…seline The five build shards each open a Gurobi session, and --processes previously defaulted to the runner's core count, opening one more session per worker process inside each shard. Both ran concurrently against the same Gurobi_Eduard WLS license, whose session baseline is 2, so the license server eventually kills sessions once the overage runs too long.
…e shards Each holds its own Gurobi session for the whole step, so running five of one concurrently already sustains more sessions than the Gurobi_Eduard WLS license's baseline of 2 for as long as the step takes; --processes further multiplied that within a single shard's scan. Splitting them into build-model and build-test lets both still run all five shards at once (no session survives from one job into the next, so the two 5-way bursts do not merge into one sustained overage), instead of serializing everything to stay under the baseline.
…ers failing needs: build-model implicitly ANDs an overall success() onto its condition, so one failed shard (build-model (4)/RPE1 hit the 90-minute job timeout in 34494319341) skipped build-test entirely, even for the four shards whose context model had already built fine.
build-model's 90-minute timeout was below what one cell line's ftINIT build can legitimately take: two extraction steps, each up to three solves (the proven-gap solve plus resolve_ties' two tie-break phases) capped at TIME_LIMIT = 30 min apiece, i.e. up to 180 min before task gap-filling. RPE1 needs about 91 min (its first step's proven-gap solve and parsimony phase both run to the 30-min cap), so its shard was cut off at the end of extraction in 34494319341 while the other four finished in 9-24 min. 240 min covers the 180-min ceiling with headroom for gap-filling and setup. The build step also went silent for the whole extraction: ftinit() only prints its step and gap-fill progress with verbose=True, and optlang sets OutputFlag = 0 on every Gurobi model it creates, so the global OutputFlag set by _set_solver_verbosity never reaches the ftINIT MILPs. verbose=True restores the step/gap-fill lines, and FTINIT_DEBUG=1 adds each solve's status and achieved gap as it finishes. Replayed locally for shard 4/5 (RPE1) with the CI's dependency versions (Python 3.11, gurobipy 13.0.3, cobra 0.32.1, raven-toolbox 61da432): prep 27 min, build-model 91 min, build-test (--processes 1) 14 min, all exiting cleanly.
The gene-essentiality workflow's aggregate job crashed with "Too many sessions, 5 active sessions for a baseline of 2": cobra prefers Gurobi automatically whenever gurobipy is importable, so loading the model for --aggregate-only opened a Gurobi WLS session even though that step never solves anything and does not need one, and the build-test matrix's 5 just-finished sessions had not yet cleared from the license server. Force glpk as the default solver before loading the model, letting only the branches that actually need Gurobi's MILP solver switch to it. Also add retry-with-backoff around that switch (a sibling shard's session is far more likely to clear within a minute or two than the license being genuinely exhausted) and explicit disposal of the default Gurobi environment at the end of each Gurobi-using branch, so a shard's session is released to the license server immediately instead of waiting on its own connection timeout.
63f1231 to
114f75a
Compare
disposeDefaultEnv() releasing a WLS session is a network round trip to the license server, unlike a local/named-user license's instant close. Running it inline risked an already-finished, already-written-to-disk shard sitting "in progress" indefinitely if that round trip stalled -- plausibly what the currently running gene-essentiality workflow run is hitting, since every shard passes through this same code regardless of how fast its own MILP solve was. Run it on a background thread with a bounded join instead: a timeout is logged and the process exits anyway, which releases the session itself, just a little later.
The gene-essentiality workflow's job logs are only retrievable after a job completes, so a currently running (or stuck) build-model/build-test shard is a black box from the outside -- no way to tell whether it is genuinely solving or wedged. Add a background resource monitor to both jobs (CPU/memory/top processes every 60s to a file, uploaded as an artifact with if: always() so it survives a timeout or failure) to give visibility into a run while or after it happens, without needing live log streaming. Also harden _use_gurobi_solver: a failed session-acquisition attempt leaves gurobipy's process-wide default environment half-initialised (the crash is inside its own construction), so best-effort dispose it before each retry rather than reusing whatever state a broken attempt left behind. Log how many attempts and how long acquisition took, plus CPU count and Gurobi's resolved Threads, separating "time spent getting a session" from "time spent solving" and surfacing possible CPU oversubscription across the 5 concurrent shards.
Main improvements in this PR:
First increment towards a pip-installable Python interface to Human-GEM (dual-language support).
human_gem.load_model()loads Human-GEM as a fully annotated cobra model — reading the YAML via raven-toolbox and merging the cross-references + SBO terms from the annotation tables, which a barecobra.io.load_yaml_modeldoes not do:pyproject.toml— packagehuman-gem(deps:cobra,pandas,raven-toolbox)human_gem/io.py—load_model(model_dir=None, annotate=True)human_gem/annotation.py— the TSV → MIRIAM/SBO merge (packaged home of the logic currently incode/annotateGEM.py)Verified locally (raven-toolbox 0.3.0):
load_model()attaches cross-references to 12698/12877 reactions, 8449/8460 metabolites and all 2848 genes;annotate=Falsereturns the bare model.This increment is deliberately small and non-disruptive — it adds new files and imports nothing existing, so it does not touch the current MATLAB or Python CI. Planned next increments:
code/annotateGEM.pyto a thin shim importing fromhuman_gem.annotation(single source of truth) + install the package in the Python CI;pip installdistribution.No model changes.
I hereby confirm that I have:
data/deprecatedIdentifiers/.developas target branch, and will be resolved with a squash-merge.mainas target branch, and will be resolved with a merge commit.