Skip to content

feat: add human_gem Python package with load_model() - #1069

Open
edkerk wants to merge 24 commits into
developfrom
feat/human-gem-python-package
Open

feat: add human_gem Python package with load_model()#1069
edkerk wants to merge 24 commits into
developfrom
feat/human-gem-python-package

Conversation

@edkerk

@edkerk edkerk commented Jul 17, 2026

Copy link
Copy Markdown
Member

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 bare cobra.io.load_yaml_model does not do:

import human_gem
model = human_gem.load_model()   # annotated cobra.Model
  • pyproject.toml — package human-gem (deps: cobra, pandas, raven-toolbox)
  • human_gem/io.pyload_model(model_dir=None, annotate=True)
  • human_gem/annotation.py — the TSV → MIRIAM/SBO merge (packaged home of the logic currently in code/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=False returns 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:

  1. reduce code/annotateGEM.py to a thin shim importing from human_gem.annotation (single source of truth) + install the package in the Python CI;
  2. promote the ftINIT gene-essentiality and metabolic-task helpers into the package;
  3. bundle/resolve the model files for a true pip install distribution.

No model changes.

I hereby confirm that I have:

  • Any removed reactions and metabolites have been moved to the corresponding deprecated identifier lists in data/deprecatedIdentifiers/.
  • This PR has develop as target branch, and will be resolved with a squash-merge.
  • This PR has main as target branch, and will be resolved with a merge commit.

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.
@github-actions

github-actions Bot commented Jul 17, 2026

Copy link
Copy Markdown

Model quality report

⚠️ 6 pre-existing finding(s), no regressions vs develop. Non-blocking.

Each check name links to its explanation in the testResults README.

Model checks

Duplicate keys (model unloadable) and no growth block the merge; every other row is a non-blocking report.

Check Result Δ vs develop
Duplicate !!omap keys 0 0
Growth (biomass producible) 125 0
Reactions with no metabolites 0 0
Model / annotation-table inconsistencies 0 0
Removed reactions or metabolites not deprecated 0 0
Metabolites missing formula 0 0
Metabolites missing charge 0 0
Reaction bound / GPR issues 0 0
Exact-duplicate reaction groups 0 0
Unused metabolites 0 0
Unused genes 0 0
Malformed cross-references 0 0
Cross-refs inconsistent across compartments 3 0 ⚠️

MACAW and mass/charge balance

Check Result Δ vs develop
Reactions flagged by MACAW dead-end test 2510 0 ⚠️
Reactions flagged as MACAW duplicates 377 0 ⚠️
Mass-imbalanced reactions 87 0 ⚠️
Charge-imbalanced reactions 234 0 ⚠️
Structure vs formula/charge inconsistencies 397 0 ⚠️

Model file and metabolic tasks

Check Result
YAML round-trip (cobrapy) pass
YAML round-trip (RAVEN) pass
YAML lint pass
Essential metabolic tasks 57 passed
Verification metabolic tasks 21 passed

MEMOTE

Total score: 63.2% (core subset)   0

Section Score Δ vs base
consistency 42.4% 0
annotation_met 73.0% 0
annotation_rxn 72.7% 0
annotation_gene 46.7% 0
annotation_sbo 81.7% 0
Per-test scores
Section Test Score
Consistency Stoichiometric Consistency 100.0%
Consistency Mass Balance 0.8%
Consistency Charge Balance 2.1%
Consistency Metabolite Connectivity 0.0%
Consistency Unbounded Flux In Default Medium 100.0%
Annotation - Metabolites Presence of Metabolite Annotation 0.0%
Annotation - Metabolites Metabolite Annotations Per Database 62.3%
Annotation - Metabolites Metabolite Annotation Conformity Per Database 45.8%
Annotation - Metabolites Uniform Metabolite Identifier Namespace 0.0%
Annotation - Reactions Presence of Reaction Annotation 0.0%
Annotation - Reactions Reaction Annotations Per Database 75.9%
Annotation - Reactions Reaction Annotation Conformity Per Database 33.3%
Annotation - Reactions Uniform Reaction Identifier Namespace 0.0%
Annotation - Genes Presence of Gene Annotation 0.0%
Annotation - Genes Gene Annotations Per Database 80.0%
Annotation - Genes Gene Annotation Conformity Per Database 80.0%
Annotation - SBO Terms Metabolite General SBO Presence 0.0%
Annotation - SBO Terms Metabolite SBO:0000247 Presence 0.1%
Annotation - SBO Terms Reaction General SBO Presence 0.0%
Annotation - SBO Terms Metabolic Reaction SBO:0000176 Presence 0.0%
Annotation - SBO Terms Transport Reaction SBO:0000185 Presence 0.7%
Annotation - SBO Terms Exchange Reaction SBO:0000627 Presence 0.0%
Annotation - SBO Terms Demand Reaction SBO:0000628 Presence 100.0%
Annotation - SBO Terms Sink Reactions SBO:0000632 Presence 100.0%
Annotation - SBO Terms Gene General SBO Presence 0.0%
Annotation - SBO Terms Gene SBO:0000243 Presence 0.0%
Annotation - SBO Terms Biomass Reactions SBO:0000629 Presence 0.0%

Full suite: 64.2%   0 · from the last /run memote.

The score above is the fast core subset. Comment /run memote to run the full suite on this pull request; the score updates here when it finishes.

Gene essentiality (Hart 2015)

Not run automatically (it takes hours). Comment /run gene-essentiality to run it on this pull request; the result posts as its own comment.

❌ = a count rose vs the target branch (regression) · ⚠️ = a pre-existing non-zero finding (non-blocking) · ⏳ = still running. Counts link to the CSV listing the exact entries.

Full workflow run · this comment is edited as results come in

@edkerk edkerk added this to the 2.1.0 milestone Jul 18, 2026
edkerk and others added 5 commits September 10, 2026 01:04
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.
The two regenerated result files conflicted; taken from develop so the
legacy getINITModel2 removal from #1068 and the graded gene-essentiality
documentation from #1077 are both kept.
@github-actions

Copy link
Copy Markdown

Gene essentiality was run on this pull request via GH Actions. Essentiality for a viability task (GR growth, ER energy/redox) is reported separately from essentiality for a capability task (SU, BS, IC), and the biomass growth ratio gives a threshold-free score (see #1076). Summary:

Graded gene essentiality vs Hart 2015 (task-scope analysis)

cellLine allTaskMCC allTaskFP viabilityMCC viabilityFP growthAUROC growthAUPRC baseRate capabilityOnly
DLD1 0.3772 127 0.3892 95 0.675 0.3382 0.1684 45
GBM 0.3208 147 0.3629 102 0.6617 0.312 0.1625 53
HCT116 0.3807 130 0.3809 103 0.6753 0.3446 0.1812 41
HELA 0.3307 164 0.3265 133 0.6696 0.2816 0.1469 44
RPE1 0.2622 174 0.3095 122 0.6585 0.2626 0.1359 59
all 0.6682 0.3081 0.1592

The full per-gene matrix is committed to data/testResults/gene-essential.csv, holding for each cell line the confusion class, the task categories the knockout breaks and the biomass growth ratio.

Note: in the case of multiple test runs, this post is edited.

…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`.
@SysBioChalmers SysBioChalmers deleted a comment from github-actions Bot Sep 10, 2026
edkerk and others added 3 commits September 10, 2026 12:55
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.
@SysBioChalmers SysBioChalmers deleted a comment from github-actions Bot Sep 10, 2026
@edkerk

edkerk commented Sep 10, 2026

Copy link
Copy Markdown
Member Author

/run gene-essentiality

@github-actions

Copy link
Copy Markdown

🧬 Started gene-essentiality on feat/human-gem-python-package (this takes a few hours). Track it in the Actions tab.

edkerk and others added 9 commits September 10, 2026 13:20
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.
@edkerk
edkerk force-pushed the feat/human-gem-python-package branch from 63f1231 to 114f75a Compare September 11, 2026 08:03
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant