Overview

This project is a small monorepo around gene set enrichment, eval exploration, and curated gene set interpretation.

The main surfaces are:

• genesets-rs: a fast Rust compute engine and CLI for ontology-aware and flat gene set enrichment.
• genesets-workflows: a Python workflow package for source preparation, configured evals, reports, metadata, curation helpers, and local tools.
• Web explorer: a local browser for generated eval report bundles.
• Curated interpretation corpus: LinkML-validated YAML interpretations of non-GO gene sets, intended both as a gold standard and as browsable content.

The repository should stay together for now. The pieces share fixtures, schemas, eval definitions, and version history, but each has a clear boundary.

Core Model

The compute engine is deliberately ontology-neutral. Gene Ontology is the main near-term target, but the core data model only requires:

• a table of term CURIEs and human-readable names;
• a precomputed child-to-ancestor closure table;
• gene-to-term annotations;
• optional gene names or symbols;
• query sets and background sets.

Flat libraries, such as MSigDB-style collections, are represented as terms with no closure. User samples are also just flat terms. This keeps the statistical engine the same for sample-vs-ontology, sample-vs-library, and arbitrary term-vs-term matrix jobs.

Boundaries

The Rust core should accept normalized tables and do one job quickly. It should not know about species, evidence codes, ontology release policy, identifier mapping, remote services, curation judgments, or web UI state.

Those concerns belong in the workflow and curation layers:

• GO, GOA, Reactome, MyGeneSet, and other source-specific adapters;
• evidence-code filters and release metadata;
• configured report runs and DuckDB summaries;
• local or static web views over generated artifacts;
• LinkML validation and curated interpretation pages.

Start with Choosing A Path if you are not sure which surface you need.

Project Map

This repository is a monorepo for four related surfaces:

The Rust crate should stay focused on the compute model. GO, GOA, Reactome, MyGeneSet, evidence-code policies, report layouts, web views, and curation workflow are deliberately outside the Rust core because they change faster and are easier to audit as files and workflow artifacts.

Repository Layout

src/
  Rust library and CLI implementation.
python/genesets-workflows/
  Python package for source prep, reports, curation helpers, and the local explorer.
evals/
  Small and parameterized eval definitions.
curation/
  LinkML schema, validator config, manifest, and curated interpretation YAML.
docs/src/
  mdBook documentation.
notebooks/
  Analysis notebooks that consume generated workflow artifacts.

Release Boundaries

Keep one repository for now, with separate release surfaces:

• genesets-rs: crates.io or cargo install release surface.
• genesets-workflows: PyPI or editable Python package release surface.
• curated interpretation corpus: versioned data release tied to repository tags.
• docs and static generated browser pages: GitHub Pages or another static host.

Split the repository only if one of those surfaces needs an independent governance, release cadence, or storage policy that makes the monorepo painful. That is not true yet.

Choosing A Path

Use this page to decide which part of the project to start with.

I Have Normalized Gene Sets And Want Enrichment

Use the Rust CLI directly:

genesets-rs matrix \
  --target-sets library.gmt \
  --target-format gmt \
  --queries samples.gmt \
  --query-format gmt \
  --output-format parquet \
  --output results.parquet

Start with Getting Started, then read the Input Model and CLI Reference.

I Need To Fetch Sources Or Run A Repeatable Report

Use genesets-workflows. The workflow layer owns downloads, source-specific filtering, GO/GOA or Reactome preparation, metadata, and report summaries:

uv run --project python/genesets-workflows genesets-workflows doctor
genesets-workflows go-impact evals/go_impact_5y_expression500.yaml

Start with Reports And Evals and Workflow Layer.

I Want To Inspect Existing Eval Results In A Browser

Use the local web explorer. It reads workflow report bundles and does not run enrichment itself:

just browser

Start with Web Explorer and Report Bundles.

I Want Curated Gene Set Interpretations

Use the curation/ corpus. Each curated YAML file records what a non-GO gene set should mean in GO terms, including curator categories, confidence, specificity, recovery status, enrichment stats, and evidence.

Start with Curated Gene Set Interpretations, then read the Curation Workflow and Schema Reference.

I Want A Public Browser For Curated Gene Sets

Build it as static pages or a workflow-layer web view over the curated YAML, not as a Rust-core feature. The recommended path is described in Static Browser Plan.

Compute Engine

The Rust compute engine is the small, fast, ontology-neutral layer. It accepts normalized tables and gene set files, performs over-representation analysis or result comparison, and writes durable output formats.

Use this layer when:

• inputs are already normalized;
• you want a single enrichment, matrix run, or result comparison;
• the task should be reproducible from explicit local files;
• workflow-specific concerns such as source download, evidence filtering, and report layout are already handled elsewhere.

The core engine does not know about species, GO release policy, evidence codes, identifier mapping, remote services, curation judgments, or web UI state. Those belong in the workflow and curation layers.

Core Pages

• Getting Started: install and run the bundled examples.
• Input Model: normalized table and gene set formats.
• CLI Reference: command surfaces and examples.
• Statistics: enrichment statistic and correction model.
• Performance Model: bitsets, matrix mode, and output costs.
• Storage Backends: TSV and Parquet output choices.

Getting Started

Build the CLI:

cargo build --release

Install this checkout so genesets-rs can be run from any directory:

cargo install --path /path/to/genesets-rs --force

Cargo writes the binary to ~/.cargo/bin/genesets-rs. If your shell cannot find it, add Cargo’s bin directory to your PATH:

echo 'export PATH="$HOME/.cargo/bin:$PATH"' >> ~/.zshrc
source ~/.zshrc

Run the bundled ontology-style example:

genesets-rs enrich \
  --annotations examples/gene_terms.tsv \
  --terms examples/terms.tsv \
  --closure examples/closure.tsv \
  --sample examples/sample.txt \
  --background examples/background.txt \
  --overlap-genes

Run all ontology terms against all ontology terms:

genesets-rs matrix \
  --annotations examples/gene_terms.tsv \
  --terms examples/terms.tsv \
  --closure examples/closure.tsv \
  --queries-from-targets \
  --background examples/background.txt

Run from YAML:

genesets-rs run examples/enrich.yaml

Relative paths inside a YAML config are resolved from the config file’s directory, so this also works from outside the repository:

genesets-rs run /path/to/genesets-rs/examples/enrich.yaml

For development without installing, use cargo run -- before the same CLI arguments:

cargo run -- run examples/enrich.yaml

Install or run the optional workflow layer when you want repeatable source prep, GO/Reactome evals, Parquet summaries, and report metadata:

uv run --project python/genesets-workflows genesets-workflows doctor

The workflow layer calls the Rust CLI for batch compute; it is not required for simple one-off enrichment.

Build the documentation site locally:

cargo install mdbook
mdbook serve

Run the test suite:

cargo test

Run the benchmark harness:

cargo bench

Input Model

The engine has one central abstraction: a named set of genes. Ontology terms, flat database entries, user samples, and backgrounds can all be represented this way.

Ontology Targets

Ontology enrichment targets are built from three tables:

Terms:

term_id    name
GO:0006915 apoptosis

Closure:

child      ancestor
GO:0006915 GO:0006915
GO:0006915 GO:0012501

Annotations:

gene_id    term_id
TP53       GO:0006915

The closure table is expected to be reflexive, but the loader defensively adds each term as its own ancestor. If a gene is annotated to a child, the prepared bitset for every ancestor receives that gene.

Flat Targets

Flat targets skip closure. They can be read from:

• list: one gene per line, one set total;
• pairwise: set_id, gene_id;
• gene-term: gene_id, set_id;
• gmt: set_id, description, genes;
• gmx: first row set IDs, remaining rows genes;
• gmx-desc: first row set IDs, second row descriptions, remaining rows genes.

Background

If no background file is supplied, the background is the union of all target genes. A supplied background file overrides that behavior and is read as one gene per row.

This matters for experimental data. For RNA-seq differential-expression lists, the background should usually be all genes that were measured and eligible to become significant, not every protein-coding gene.

Gene Names

Gene names are optional. The statistical engine uses gene IDs only. A gene-name table is only used for human-readable overlap output.

CLI Reference

The CLI has four subcommands:

• enrich: one query set against many target sets;
• matrix: many query sets against many target sets;
• run: YAML-configured enrich or matrix.
• compare: threshold-crossing diff between two result tables.

enrich

genesets-rs enrich \
  --annotations gene_terms.tsv \
  --terms terms.tsv \
  --closure closure.tsv \
  --sample sample.txt \
  --background background.txt \
  --output results.tsv

Use --target-sets instead of --annotations for flat libraries:

genesets-rs enrich \
  --target-sets library.gmt \
  --target-format gmt \
  --sample sample.txt

Useful options:

• --sample-format: auto, list, pairwise, gene-term, gmt, gmx, gmx-desc;
• --sample-set: select one set from a multi-set sample file;
• --min-overlap: suppress rows below an overlap count;
• --max-p-value: suppress rows above a raw p-value cutoff;
• --max-p-adjust: suppress rows above an adjusted p-value cutoff;
• --correction: bonferroni or none;
• --output-format: tsv, parquet, or null; parquet requires --output, and null is useful for compute-only profiling;
• --overlap-genes: include overlapping gene IDs and names;
• --threads: set Rayon worker count.

matrix

genesets-rs matrix \
  --annotations gene_terms.tsv \
  --terms terms.tsv \
  --closure closure.tsv \
  --queries queries.gmx \
  --query-format gmx \
  --background background.txt

Use targets as queries for term-vs-term runs:

genesets-rs matrix \
  --annotations gene_terms.tsv \
  --closure closure.tsv \
  --queries-from-targets

run

mode: enrich
ontology:
  terms: terms.tsv
  closure: closure.tsv
  annotations: gene_terms.tsv
input:
  sample: sample.txt
  sample_format: list
  sample_name: sample
background:
  file: background.txt
overlap_genes: true
max_p_adjust: 0.05
output_format: tsv

Run it:

genesets-rs run examples/enrich.yaml

Relative paths in YAML are resolved from the config file’s directory. Relative paths passed directly as CLI arguments are resolved from the current working directory.

For mass evals, write Parquet and inspect it with DuckDB:

genesets-rs matrix ... --output-format parquet --output results.parquet
duckdb -c "SELECT * FROM 'results.parquet' WHERE p_adjust_bonferroni <= 0.05"

compare

Compare two enrichment result tables by (query_id, target_id) and classify adjusted p-value threshold crossings:

genesets-rs compare \
  --left go-2021.parquet \
  --right go-2026.parquet \
  --p-adjust-cutoff 0.05 \
  --output-format parquet \
  --output go-2021-vs-2026.diff.parquet \
  --metadata-output go-2021-vs-2026.diff.yaml

Input formats are inferred from .tsv, .txt, .parquet, or .pq, or can be set explicitly with --left-format and --right-format. TSV output goes to stdout by default. Parquet output requires --output.

The default output includes:

• lost_significant: significant on the left, not significant on the right;
• gained_significant: not significant on the left, significant on the right;
• shared_significant: significant on both sides.

Use --crossings-only to emit only gained/lost rows.

Workflow CLI

The separate genesets-workflows command is the convenience layer for configured source prep and reports. It calls the Rust CLI for batch compute and then writes Parquet plus metadata:

uv run --project python/genesets-workflows genesets-workflows doctor
genesets-workflows go-impact evals/go_impact_5y_expression500.yaml
genesets-workflows reactome-flat

go-impact YAML can apply query_sets.include_regex and query_sets.exclude_regex before selecting query sets from a cached snapshot. The regexes match query id, name, description, source class, stratum, and search query metadata. For example, exclude GO-derived MSigDB query sets from a GO target analysis:

query_sets:
  source_dir: expression_like/generated/msigdb_diverse_5k
  limit: 4313
  exclude_regex:
    - "^(GOBP|GOCC|GOMF)_"

For interactive review of an existing report bundle, use the optional web explorer dependencies:

uv run --project python/genesets-workflows --extra explorer \
  genesets-workflows explore notebooks/generated/go_iba_impact_expression5000_diverse

Use genesets-rs for normalized single jobs. Use genesets-workflows when the task needs downloads, evidence filters, release metadata, multiple Rust runs, DuckDB summaries, notebook/report artifacts, or browser-based result triage.

Statistics

The MVP performs standard one-sided over-representation analysis.

For each query and target term, the engine computes:

The p-value is the right-tail hypergeometric probability:

P(X >= k), X ~ Hypergeometric(N, K, n)

This is equivalent to one-sided Fisher exact enrichment on the 2 x 2 table:

Multiple Testing

The current correction choices are:

• bonferroni: min(p * number_of_tests, 1);
• none: report raw p-values as adjusted p-values.

For matrix runs, the Bonferroni denominator is the number of non-empty query by non-empty target tests.

Future Statistics

The code is structured so additional test families can be added without changing the data model:

• Benjamini-Hochberg FDR;
• ranked-list methods;
• topology-aware parent-child methods;
• elim and weight-style algorithms;
• model-set approaches inspired by Ontologizer.

Performance Model

genesets-rs trades initialization work and memory for fast repeated calculations.

Bitset Representation

Each query and target is represented as a dense Vec<u64> over a shared gene universe. After loading:

• set size is a popcount;
• overlap is a word-wise AND plus popcount;
• ontology propagation is already baked into target bitsets;
• matrix mode parallelizes over query sets.

This is a good fit when the universe is moderate, memory is available, and the same target library is reused across many query sets.

Expected Bottlenecks

There are three separate costs to keep apart:

• preparing inputs: parsing term, closure, annotation, background, and query files, then building dense bitsets;
• scoring: intersecting query and target bitsets, running Fisher exact tests, and applying cutoffs;
• writing results: resolving string labels and serializing rows.

The CLI has --output-format null for compute-only timing. Use it to separate scoring from serialization:

genesets-rs matrix ... --output-format null

On the current 200 MSigDB GSE-style query fixture against GOA human and the current GO closure:

• one query, significant-only, no output serialization: about 1.4s with one Rayon thread;
• 200 queries, significant-only, no output serialization: about 5.9s with one Rayon thread, or about 2.2s on a 10-core laptop;
• 200 queries, no p-value cutoff, no output serialization: about 2.25s wall time and about 540 MB peak RSS;
• 200 queries, no p-value cutoff, TSV output: about 3.6s wall time for about 778k result rows and a 140 MB TSV.

For significant-only evals, TSV is not the main bottleneck. Rebuilding the ontology target index on every CLI process dominates wall time. For unfiltered or weakly filtered mass comparisons, row materialization and TSV serialization become large enough to matter.

TSV remains useful for individual analyses, but multi-timepoint and many-query evals should move to a columnar or embedded analytical format. The output module is intentionally isolated from the scoring layer so Parquet and DuckDB sinks can be added without changing enrichment statistics.

Current Diagnosis

The most important bottlenecks are:

• repeated text parsing and ontology target bitset construction for every run;
• duplicated string-heavy closure and annotation tables during setup;
• materializing every retained result row in memory before writing;
• sorting all retained rows before output;
• dense bitset scans over every query-target pair, even for very sparse terms;
• TSV formatting when result sets are large or unfiltered.

The next performance milestones are:

• serialized prepared indexes for term bitsets, term metadata, and gene metadata;
• streaming or batched result sinks so large runs do not materialize every row;
• Parquet output for columnar scans and archival eval artifacts;
• DuckDB output for direct diff queries over many timepoints;
• top-k and p-value-threshold output modes;
• mmap-friendly bitset storage;
• sparse or Roaring-style bitsets benchmarked against dense bitsets on real GO fixtures;
• criterion benchmarks on synthetic and real ontology fixtures;
• comparison runs against GOATOOLS, g:Profiler, topGO, and clusterProfiler on the same frozen inputs.

Result Storage

For the mass-diff use case, the preferred storage shape is numeric and columnar:

• dimension tables for runs, queries, targets, terms, genes, ontology releases, and annotation variants;
• fact rows keyed by (run_id, query_index, target_index);
• numeric columns for overlap, sizes, p-value, adjusted p-value, and significance flags;
• optional overlap-gene payloads in a separate table or sidecar file.

Parquet is a good artifact format for immutable eval output and batch scans. DuckDB is a good working format for interactive diff queries, threshold-crossing queries, joins against ontology metadata, and summaries across many releases. Both should sit behind the output boundary rather than leak into the enrichment model.

Benchmark Harness

The repository includes a Criterion benchmark:

cargo bench

The benchmark is synthetic by design. It measures the core scoring path without network access, identifier mapping, or ontology parsing.

Storage Backends

The result backend should optimize for fast writes, compact artifacts, and easy diff queries over many runs. TSV should remain available for small, inspectable single analyses, but it should not be the mass-eval format.

Recommendation

Use Parquet as the first durable mass-eval output format. Use DuckDB as the default analysis layer over Parquet files rather than making .duckdb the first primary output format.

This gives us:

• one fast, compact artifact format for CI, eval archives, and sharing;
• direct SQL introspection through DuckDB with no import step;
• simpler Rust output code using Arrow record batches and the native Parquet writer;
• a clean path to add a DuckDB sink later if repeated interactive sessions need a materialized database.

Supporting both as eventual output options is reasonable, but implementing both as first-class write paths immediately is likely overkill. The first target should be:

genesets-rs ... --output-format parquet --output results.parquet
duckdb -c "SELECT * FROM 'results.parquet' WHERE p_adjust_bonferroni <= 0.05"

Diff outputs use the same approach:

genesets-rs compare --left old.parquet --right new.parquet \
  --output-format parquet --output old-vs-new.diff.parquet
duckdb -c "SELECT class, count(*) FROM 'old-vs-new.diff.parquet' GROUP BY class"

Comparative Notes

DuckDB can read Parquet directly and can push projections and filters into the Parquet scan. That means we do not need to choose between Parquet artifacts and DuckDB introspection.

Local Smoke Test

Using the 200 MSigDB GSE-style query fixture against current GOA human, the unfiltered TSV had 778,450 result rows.

DuckDB query timings, including process startup, were effectively tied on this fixture:

These numbers are too small to settle large-scale behavior, but they are enough to reject the idea that Parquet would cost us DuckDB-style introspection.

The first Rust Parquet backend is optimized for compatibility and write speed with Snappy compression. On the same fixture, end-to-end unfiltered output was about 3.7s for Parquet versus about 3.5s for TSV, while reducing output size from 132 MB to 11 MB. The synthetic writer benchmark is more favorable to Parquet: about 1.45 ms for Parquet versus about 3.7 ms for TSV on the current 150 x 750 fixture.

Implementation Direction

The output boundary should produce numeric, batched result records:

• run_id, query_index, target_index;
• overlap, query_size, target_size, background_size;
• p_value, p_adjust, and boolean significance flags;
• optional string dictionaries or dimension tables for query and target labels;
• optional overlap genes in a separate sidecar table/file.

The Parquet backend writes Snappy-compressed row groups from these batches. A later DuckDB backend can consume the same batches through an appender or materialize a DuckDB database from the Parquet outputs.

Reports And Evals

Reports and evals live above the Rust compute engine. They answer questions that require source-specific preparation, repeated runs, metadata, comparison logic, and human review.

The main implementation lives in python/genesets-workflows. It prepares inputs, calls genesets-rs in coarse batches, writes Parquet and YAML/JSON artifacts, and then summarizes those artifacts with DuckDB, notebooks, docs, or the local web explorer.

What Belongs Here

• fetching public source data from GO, GOA, Reactome, MyGeneSet, or similar;
• applying evidence-code and relation policies;
• preparing normalized term, closure, annotation, query, and background files;
• running configured batches through the Rust CLI;
• comparing ontology releases or annotation variants;
• writing report bundles that can be inspected by SQL, notebooks, or the web UI;
• using curated interpretations as precision/recall fixtures.

Core Pages

• Workflow Layer: package boundary and commands.
• Evals: current eval datasets and report findings.
• Diffing: threshold-crossing comparisons.
• Post-Processing: term filtering and interpretation.
• Tutorials: runnable shell and notebook workflows.
• Notebook Workflows: notebook conventions.

Workflow Layer

genesets-rs is the compute kernel. It should stay small, ontology-neutral, and fast. The outer convenience layer lives in the Python package under python/genesets-workflows.

The workflow layer owns the parts that are useful but not core enrichment logic:

• fetching public sources such as GO, GOA, Reactome, and MyGeneSet;
• preparing normalized tables, closures, evidence-code variants, and backgrounds;
• running batch analyses through the Rust CLI;
• writing Parquet result artifacts and YAML metadata;
• querying result Parquet with DuckDB;
• computing global report diagnostics such as GO term coverage and terms that are scorable but never appear as significant enrichment hits;
• generating report summaries for docs, notebooks, and web views;
• validating curated interpretations and, eventually, generating curated corpus browser pages.

Local Use

Run the packaged workflow CLI from this checkout:

uv run --project python/genesets-workflows genesets-workflows doctor

For repeated local use, install it editable:

python3 -m pip install -e python/genesets-workflows
genesets-workflows doctor

The Rust binary still needs to be installed or on PATH:

cargo install --path . --force

From outside the repository, use an absolute project path or install the workflow package into the active environment:

uv run --project /path/to/genesets-rs/python/genesets-workflows \
  genesets-workflows doctor

python3 -m pip install -e /path/to/genesets-rs/python/genesets-workflows

When To Use Which Layer

Use genesets-rs directly when inputs are already normalized and the task is a single enrichment, matrix enrichment, or result comparison:

genesets-rs matrix \
  --target-sets reactome.gmt \
  --target-format gmt \
  --queries expression_sets.gmt \
  --query-format gmt \
  --output-format parquet \
  --output results.parquet

Use genesets-workflows when the task has source-specific prep, multiple Rust runs, metadata, or report summaries:

genesets-workflows go-impact evals/go_impact_5y_expression500.yaml

The workflow command should usually call Rust once per large batch, not once per gene set. If a workflow needs thousands of enrichments, it should express that as a batch plan for the Rust CLI.

Launch the local web explorer over an existing report bundle when the task is interactive triage rather than batch computation:

just browser

The default browser recipe loads the current GOA all-vs-IBA report, the 2021-vs-2026 GO/GOA temporal report, and the current GOA all-vs-no-contributes_to report. Use just browser-iba, just browser-go5y, or just browser-contributes to open one analysis directly.

uv run --project python/genesets-workflows --extra explorer \
  genesets-workflows explore notebooks/generated/go_iba_impact_expression5000_diverse

Use the curation commands when the task is adjudicating or validating curated GO interpretations:

uv run --project python/genesets-workflows --extra curation \
  genesets-workflows curate --help

Global Diagnostics

The go-impact report also writes term-coverage Parquet files. These join the prepared GO term table, propagated annotation target sizes, and retained enrichment results. Each term is classified as:

• significant: appears in at least one retained enrichment row;
• scorable_never_significant: has propagated annotations in the analysis background but never crosses the report threshold for the query collection;
• unscorable: has no propagated annotations in the analysis background.

This answers a different question from threshold-crossing diffs: which parts of the ontology are effectively invisible for this query collection and annotation variant? The paired coverage file compares those statuses across the two snapshots or annotation variants.

Current Commands

GO impact report over two prepared GO snapshots:

genesets-workflows go-impact evals/go_impact_5y_expression500.yaml

Expression signatures against official Reactome as a flat pathway library:

genesets-workflows reactome-flat

Prepare the official Reactome GMT target library:

genesets-workflows prepare-reactome-flat \
  --out-dir evals/reactome_flat/generated/current

Fetch a MyGeneSet query snapshot:

genesets-workflows fetch-mygeneset \
  --query 'GSE*' \
  --source-filter msigdb \
  --limit 500 \
  --out-dir evals/expression_like/generated/msigdb_gse_500

For report-quality eval fixtures, prefer a stratified source plan over the first GSE* hits:

genesets-workflows fetch-mygeneset-stratified \
  evals/expression_like/msigdb_diverse_5k.yaml

The stratified fetcher still emits ordinary queries.gmt, background.txt, and metadata.json files. The difference is that metadata.json records the quota stratum and search query for each set, which makes later report examples auditable.

go-impact configs can also filter the selected query snapshot with regexes before running Rust. This is useful for benchmark hygiene, for example excluding GO-derived MSigDB query sets when the target ontology is GO:

query_sets:
  source_dir: expression_like/generated/msigdb_diverse_5k
  limit: 4313
  exclude_regex:
    - "^(GOBP|GOCC|GOMF)_"

The report metadata records the selected, considered, and skipped query counts.

Browse one or more generated report bundles in a local web UI:

genesets-workflows explore notebooks/generated/go_iba_impact_expression5000_diverse

The old script entry points for these packaged commands remain as compatibility wrappers around the package modules. Some older eval helpers are still script-native and should be migrated as the workflow layer expands.

Interface Boundary

Python should call Rust through batch-oriented commands, not tight subprocess loops. The desired pattern is:

Python source prep -> normalized tables/YAML plan -> one Rust batch command
-> Parquet outputs -> DuckDB summaries/reports

This keeps the subprocess boundary coarse. If a future notebook, web service, or Python API needs in-process enrichment, a maturin binding can wrap the same Rust core crate without replacing the CLI.

The current choice is CLI-first rather than maturin-first because the most variable code is not the Fisher test. It is fetching, filtering, release metadata, identifier normalization, report generation, and SQL over Parquet. Those pieces are faster to evolve in Python. The Rust interface should grow batch abstractions before Python grows tight bindings.

Artifact Layout

Workflow outputs should be predictable so notebooks, docs, CI jobs, and web viewers can consume the same products:

run-dir/
  summary.yaml
  summary.json
  queries.gmt
  queries.metadata.json
  left-results.parquet
  right-results.parquet
  left-vs-right.diff.parquet
  left-vs-right.diff.yaml

Use Parquet for durable result tables. Use small YAML/JSON files for metadata, parameters, source URLs, file hashes, version labels, row counts, and timing. Use TSV only for small human-inspectable runs.

DuckDB is the default introspection layer over Parquet:

duckdb -c "SELECT class, count(*) FROM 'left-vs-right.diff.parquet' GROUP BY class"

That keeps artifacts immutable and easy to share while preserving SQL inspection.

Distribution Plan

The repository should remain a monorepo:

• Rust crate and CLI: crates.io release surface.
• genesets-workflows: PyPI release surface.
• Docs, eval configs, and notebooks: same git history and version tags.

The Python package should record the expected genesets-rs CLI version and check it at runtime as the workflow layer matures. Release tags can drive both crates.io and PyPI publishing.

Short-term distribution:

cargo install --path .
uv run --project python/genesets-workflows genesets-workflows doctor

Longer-term distribution:

cargo install genesets-rs
python3 -m pip install genesets-workflows

If in-process Python calls become necessary, add maturin bindings around the same Rust core crate. The CLI should remain supported because it is the most transparent interface for batch runs, notebooks, workflow engines, and remote execution.

Notebooks And Web Explorer

Notebooks should demonstrate the CLI and analyze generated artifacts. They should not be the primary workflow engine. A notebook should run a configured workflow command, then display Parquet-derived tables and plots.

The web explorer follows the same model: select a run directory, read summary.yaml, query Parquet with DuckDB, and render threshold crossings, timing, top changed targets, enrichment rows, and query genes. No web-specific behavior enters the Rust enrichment kernel.

The curated gene set browser should follow the same separation. Validate the YAML corpus, materialize a JSON read model, and generate static pages or a client-side browser from that model. Do not move curation or browser behavior into the Rust enrichment core.

Evals

The eval framework should answer two questions:

1. Are results statistically consistent with trusted implementations on frozen inputs?
2. Is the engine faster on workloads we care about?

Layers

Unit tests validate local pieces: bitsets, Fisher exact p-values, loaders, and annotation propagation.

CLI integration tests validate the executable surface and TSV output shape.

Eval datasets validate end-to-end workflows. These should be explicit, versioned, and reproducible.

Benchmarks validate speed on synthetic and real workloads.

Expression20 Starter Eval

The preferred real-data smoke panel is evals/expression20. It names 20 human MSigDB/MyGeneset.info expression-derived signatures selected across:

• in vivo or primary immune cell-type contrasts;
• patient or clinical disease-state contrasts;
• vaccination and infection-response contrasts;
• cytokine-stimulated primary cells;
• drug response and resistance signatures;
• cell-line stimulation/differentiation signatures.

Fetch the current GMT snapshot:

python3 scripts/fetch_mygeneset_eval.py \
  --manifest evals/expression20/sets.tsv \
  --out-dir evals/expression20/generated

Run expression-vs-expression overlap enrichment:

genesets-rs run evals/expression20/config.yaml

Run Expression20 against official GO:

python3 scripts/run_disease20_go_eval.py \
  --manifest evals/expression20/sets.tsv \
  --out-dir evals/expression20_vs_go/generated \
  --eval-name expression20_vs_go \
  --description "Twenty expression-derived human MSigDB/MyGeneset.info signatures enriched against official GOA human GO annotations."

Disease20 Legacy Eval

Disease20 is retained as a legacy/artificial knowledge-curation contrast. It names 20 human Disease Ontology gene sets from MyGeneset.info. The manifest is intentionally small and hand-reviewed, but it is not representative of differential-expression query workflows.

Fetch the current GMT snapshot:

python3 scripts/fetch_mygeneset_eval.py \
  --manifest evals/disease20/sets.tsv \
  --out-dir evals/disease20/generated

This writes:

• queries.gmt: disease gene sets in GMT format;
• metadata.json: source URL, fetch time, gene counts, and skipped IDs.

Generated files are not required for cargo test. The goal is to keep normal CI deterministic while making real-data evals one command away.

Disease20 vs GO

The second eval uses the same 20 MyGeneset.info disease query sets against official GOA human GO annotations.

Run it:

python3 scripts/run_disease20_go_eval.py

This downloads official GO inputs, prepares ontology-neutral tables, runs four annotation-filter variants, and writes companion metadata:

• all: all evidence codes, NOT filtered, contributes_to retained;
• no_contributes_to: all evidence codes, NOT and contributes_to filtered;
• iba: IBA only, NOT filtered;
• iba_iea: IBA or IEA only, NOT filtered.

The prep stage writes metadata.yaml with source URLs, file digests, GAF header details, closure relation policy, and annotation filter counts. The run stage writes run_metadata.yaml with command timings and result summaries.

These filters are intentionally outside the Rust core engine. The core consumes prepared gene_id, term_id tables; eval helpers decide how to derive those tables from GAF.

By default, generated eval configs include max_p_adjust: 0.05, so result TSVs contain Bonferroni-significant rows. The companion metadata records this cutoff and the per-variant CLI runtime.

GO Timepoints

The go_timepoints eval runs Disease20 against official GO archive releases:

• 2026-03-25;
• 2021-05-01;
• 2016-05-01.

Run it:

python3 scripts/run_go_timepoints_eval.py

The summary metadata contains per-timepoint ontology sizes, GAF metadata, per-variant row counts, run timing, and pairwise diffs over (query_id, target_id) result keys.

Expression500 GO All vs IBA

The IBA-impact eval compares current GOA human all-evidence annotations against the IBA-only annotation variant over the same 500 expression-derived query sets used by the GO impact report:

genesets-workflows go-impact evals/go_iba_impact_expression500.yaml

This uses the same prepared GO ontology, closure, background, and query sets on both sides. Only the annotation table changes:

• left side: all, NOT-filtered GOA human, all evidence codes, contributes_to retained;
• right side: iba, NOT-filtered GOA human, IBA evidence only.

The local June 2026 run used the current 2026-03-25 GO/GOA prep. The annotation variant sizes were:

At Bonferroni-adjusted p <= 0.05, the report produced:

The threshold-crossing diff classified:

The same run now records GO term coverage, which asks which propagated GO targets are scorable but never become significant for this query collection:

The paired coverage comparison classified 10,338 terms as scorable but never significant in both variants, 9,536 as scorable but never significant on only one side, and 1,176 as significant in at least one variant. This is useful ontology “dark matter” for a fixed query collection: terms can be fully annotated and still never show up in enrichment output.

This supports the expected contraction from all evidence to IBA-only, but it also shows why post-processing matters. Because IBA changes target sizes and the correction universe, some very broad terms become significant only in the IBA run. With a simple specificity guard of target size <= 1,000, the diff becomes 3,699 lost, 1,401 retained, and 350 gained rows.

Interpretation: IBA is a useful vetted subset, but not an automatic “less-biased answer set.” It is also a lower-coverage subset, and it can shift enrichment toward terms whose IBA target sizes shrink dramatically. Reports should therefore present all-vs-IBA threshold crossings together with target size deltas and a configurable slim/antislim or size-based term filter.

Qualitative spot checks support that caution:

The working interpretation is therefore: IBA is more curated and useful as a contrast set, but many “lost” calls are plausible biology that lacks IBA coverage. Many “gained” calls are not discoveries; they are threshold crossings caused by smaller IBA target sets and a changed correction universe. The report should expose target-size deltas and representative overlap genes before promoting a gained/lost term as biologically meaningful.

The losses are not all equivalent. Useful report examples should separate desirable losses from concerning losses:

Good losses, where IBA removes broad or weakly interpretable all-evidence calls:

Bad losses, where IBA removes biologically plausible calls for the query set:

The report should therefore label losses as “possibly useful pruning” only when they are broad, recurrent, and replaceable by more specific retained terms. Losses that match the experimental contrast and have coherent overlap genes should be highlighted as IBA coverage gaps.

Some bad losses have partial compensation, but the compensation is usually from retained IBA-significant neighbors rather than from newly gained IBA-only terms:

These spot checks are still too immune-heavy. That is partly a fixture problem: the first 500 MyGeneset/MSigDB GSE* hits are dominated by C7-style immunologic signatures. The IBA report should therefore avoid presenting those examples as representative of all expression data. Less immune-centric checks from the same run are more mixed:

This suggests a useful display rule: for each lost term, show whether the same query still has retained or gained terms in the same ontology neighborhood. A loss with strong neighboring retained terms is less damaging than a loss where the entire biological neighborhood disappears.

For the next expression-scale fixture, sample by source family or keywords rather than taking the first GSE* results. The eval set should deliberately include immune, cancer, stromal, epithelial, neuronal, developmental, cell-cycle, drug-response, and stress-response signatures so qualitative examples are not all drawn from immune contrasts.

The stratified 5,000-set plan lives at evals/expression_like/msigdb_diverse_5k.yaml and can be fetched with:

genesets-workflows fetch-mygeneset-stratified \
  evals/expression_like/msigdb_diverse_5k.yaml

The same source can drive temporal GO and all-vs-IBA reports:

genesets-workflows go-impact evals/go_impact_5y_expression5000_diverse.yaml
genesets-workflows go-impact evals/go_iba_impact_expression5000_diverse.yaml

Benchmark composition matters. The refreshed source reserves modern MSigDB families rather than relying on first-N GSE* or old pathway hits:

This is a benchmark-suite rule, not just a convenience for one notebook: large gene-set benchmarks should report source-family composition and reserve space for newer collections such as MSigDB C8, C9, and 3CA.

A local synthetic scale probe used the current 500-set fixture repeated with unique query IDs. This is not a biological eval, but it answers the engineering question about 10x query-set scale:

The full 5,000-set synthetic all-vs-IBA path was about 31s on this laptop. The all-evidence side wrote 208,220 significant rows; the IBA side wrote 99,290; the diff wrote 237,560 rows. The row count is not exactly 10x the 500-set run because Bonferroni becomes stricter as query count increases.

The actual stratified 5,000-set all-vs-IBA report completed in 52.6s after refreshing the query source to include explicit C8/C9/3CA coverage:

The diff classes were 337,357 lost, 136,246 shared, and 27,620 gained significant pairs. This is a better scale signal than the repeated synthetic run because the query collection is actually diverse and has a broader gene background.

The same stratified source against the five-year temporal GO/GOA diff completed in 47.3s:

The temporal diff classes were 73,369 lost, 215,638 shared, and 257,965 gained significant pairs. This is the current best estimate for a coarse two-snapshot GO impact report over 5,000 diversified expression-like query sets.

Expression-Like Query Sets

Disease Ontology gene sets are convenient but artificial. For scale and realism, use MyGeneset/MSigDB signatures:

genesets-workflows fetch-mygeneset \
  --query 'GSE*' \
  --source-filter msigdb \
  --limit 2000 \
  --out-dir evals/expression_like/generated/msigdb_gse_2k

MSigDB GSE, *_UP, *_DN, and vs signatures are closer to differential-expression outputs than disease ontology gene sets. These should become the main large-query workloads for timing and GO-version diffing.

Expression500 vs Reactome Flat

Reactome can be tested without any new Rust code because the core engine already supports flat GMT target libraries. The current flat Reactome eval uses the same 500 expression-derived query sets as the GO impact report and enriches them against the official Reactome pathway GMT.

Run it:

genesets-workflows reactome-flat

This performs two steps:

1. genesets-workflows prepare-reactome-flat downloads ReactomePathways.gmt.zip from the official Reactome download directory and normalizes it from name, id, genes... to id, name, genes....
2. genesets-rs run evals/reactome_flat/config.yaml enriches the expression query sets against that flat pathway library.

The local May 2026 run prepared 2,557 Reactome pathways over 11,963 gene symbols and produced 509 Bonferroni-significant rows across 85 query sets. The enrichment step took about 0.22 seconds on this laptop.

This is intentionally the flat baseline. A hierarchy-aware Reactome workflow should prepare the standard ontology-style trio:

• terms.tsv from ReactomePathways.txt;
• closure.tsv from ReactomePathwaysRelation.txt;
• gene_terms.tsv from a lowest-level mapping such as NCBI2Reactome.txt, with a gene ID strategy matched to the query sets.

Planned Comparison Evals

The most useful comparison suite will freeze:

• ontology release;
• association release;
• closure relation policy;
• gene universe;
• query gene sets;
• exact correction method;
• expected significant terms from each comparator.

Comparator adapters should run each external tool in containers where possible and emit normalized TSV for comparison.

Diffing

Diffing should operate on canonical, unpostprocessed or minimally filtered result tables. Post-processing can make a term disappear for reasons unrelated to enrichment change, which is useful for presentation but bad for version comparison.

Use Cases

Common comparisons:

• GO release A vs GO release B;
• contributes_to retained vs filtered;
• IBA-only vs all evidence;
• IBA plus IEA vs all evidence;
• alternative backgrounds;
• alternative correction methods.

Stable Join Keys

The default diff key should be:

query_id, target_id

For ontology version comparisons, classify missing terms explicitly:

• present_both;
• missing_left;
• missing_right;
• obsolete_or_removed;
• merged_or_replaced, when ontology metadata can prove it;
• filtered_by_p_cutoff, when a term exists in the full run but not in the cutoff TSV.

This means serious diffs need either full result tables or at least a top-k/threshold run that is generous enough not to hide important deltas.

Delta Metrics

Useful columns:

Rank deltas are often more interpretable than tiny p-value deltas at extreme significance.

For thresholded result TSVs, the primary endpoint should be significance crossing:

A term moving from 1e-10 to 1e-20 is usually less important than a term crossing the cutoff. Large log-p deltas are still useful as a secondary diagnostic, especially when target size or background size changed dramatically.

Ordering

Recommended pipeline:

1. Prepare each ontology/annotation variant.
2. Run enrichment.
3. Diff canonical results.
4. Optionally post-process each side for display.
5. Optionally diff reduced display tables, but label that as a presentation diff.

This ordering keeps version comparisons honest. A term that disappears due to redundancy pruning should not be confused with a term that disappeared because the ontology or annotations changed.

CLI Workflow

genesets-rs compare implements the coarse threshold-crossing diff. It reads TSV or Parquet result tables, joins on (query_id, target_id), and emits rows where either side is significant:

genesets-rs compare \
  --left go-2021.parquet \
  --right go-2026.parquet \
  --p-adjust-cutoff 0.05 \
  --output-format parquet \
  --output go-2021-vs-2026.diff.parquet \
  --metadata-output go-2021-vs-2026.diff.yaml

Missing rows are treated as not significant. This makes the command usable with both full matrices and result tables that were already filtered by --max-p-adjust. If the inputs are filtered, left_present=false or right_present=false means “not present in this result table”; it does not by itself prove that the term disappeared from the ontology.

For a strict crossing-only report:

genesets-rs compare \
  --left go-2021.parquet \
  --right go-2026.parquet \
  --crossings-only \
  --output-format parquet \
  --output crossings.parquet

Metadata

Every diff output should include companion metadata:

• left/right result file paths and digests;
• left/right ontology and annotation metadata paths;
• cutoff policy;
• whether missing rows were looked up in a full table;
• diff metrics and ranking policy;
• post-processing status.

The current --metadata-output YAML includes input paths, inferred formats, cutoff policy, input row counts, output row count, and class counts. It does not yet include file digests or ontology metadata joins.

Post-Processing

GO enrichment routinely returns many related terms. The scoring engine should stay simple and complete; reduction should be a composable post-processing layer that consumes enrichment TSV plus closure/annotation metadata.

Standard Families

Significance cutoff is the first reduction step. PANTHER-style APIs expose a cutoff parameter; when FDR correction is used, the cutoff is applied to FDR, otherwise to p-value. Our eval configs mirror this with max_p_adjust: 0.05.

GO slims are a curated reduction strategy. Instead of pruning enriched terms after scoring, project annotations or results onto a smaller subset of broad terms. This is useful for overviews, less useful when the user wants precise mechanistic terms.

Slim-descendant filters are a related report-ranking strategy. Instead of projecting results onto slim terms, keep candidate terms that are themselves in a slim or descend from a slim term, often with a target-size guard. This is useful for “largest changed term” summaries where root-like terms such as broad molecular function categories would otherwise dominate.

Antislim filters invert that idea: exclude a slim region and its descendants from a report. This is useful when a known broad ontology branch is technically correct but distracts from the biological question being reviewed.

Ancestor-descendant pruning uses the closure graph. A simple policy is: for each query, if a term has a significant descendant with equal or better adjusted p-value, hide the ancestor. This is easy to explain, but it can over-prune when a broad parent captures a coherent process and many children are sparse.

Gene-overlap pruning compares term gene sets. If two enriched terms share nearly the same query-overlap genes, keep one representative. Representative selection can rank by adjusted p-value, specificity, overlap size, or information content.

GO Trimming-style approaches remove redundant parent terms based on graph paths and result statistics after enrichment has already been run.

REVIGO-style approaches cluster terms by semantic similarity and keep representatives. This is more flexible than strict ancestor-descendant pruning because related terms need not sit on the same direct path.

Topology-aware scoring methods such as topGO elim, weight, weight01, and parentchild, or Ontologizer-style parent-child/model-set methods, are not just post-processing. They change the p-values by accounting for the GO graph during scoring. We should support them later as alternate statistical models, not as TSV filters.

Proposed Composable Reducers

Reducers should take a full result table and emit:

• a reduced table;
• a mapping from hidden term to representative term;
• a reason code;
• parameters and input file digests in metadata.

Candidate reducer modes:

The tolerance is important. Sometimes a general parent has a slightly better p-value because it adds many relevant genes. A practical rule is not binary dominance, but dominance with a margin:

hide ancestor A if descendant D is significant and
  p_adjust(D) <= p_adjust(A) * tolerance

With tolerance = 1.0, the descendant must be at least as significant. With tolerance = 2.0, a slightly worse but much more specific descendant can still represent the signal.

PANTHER Behavior

PANTHER’s public service exposes Fisher or binomial tests, FDR/Bonferroni/none correction, and a cutoff parameter. Public documentation emphasizes cutoff-filtered result tables and GO-slim datasets. I did not find evidence that PANTHER applies a dynamic ancestor-descendant redundancy-pruning pass to complete GO over-representation outputs. For our purposes, treat PANTHER as a reference for cutoff filtering and GO-slim-style summaries, not as a specific redundancy-pruning algorithm.

Recommended Default

For eval tables, keep the unpruned significant TSV as the canonical result:

max_p_adjust: 0.05
post_processing: none

Then create optional reduced views:

post_processing:
  method: ancestor-dominated
  p_adjust_tolerance: 2.0
  prefer: specificity

This avoids losing terms before diffing or debugging.

Sources

• PANTHER API parameters via rbioapi docs: https://rbioapi.moosa-r.com/reference/rba_panther_enrich.html
• GO subset guide: https://geneontology.org/docs/go-subset-guide/
• GO Trimming paper: https://pmc.ncbi.nlm.nih.gov/articles/PMC3160396/
• REVIGO paper: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0021800
• topGO manual: https://bioconductor.org/packages/devel/bioc/vignettes/topGO/inst/doc/topGO_manual.html

Tutorials

These tutorials cover compute and report workflows. Each one has a matching Jupyter notebook under the repository’s notebooks/ directory, but the commands shown here can also be run directly in a shell.

Notebook Setup

Install the CLI from the checkout:

cargo install --path /path/to/genesets-rs --force

Install the workflow/notebook environment from the Python package:

uv run --project python/genesets-workflows --extra notebooks genesets-workflows doctor

Generated notebook outputs should be written under notebooks/generated/, which is ignored by git.

Available Tutorials

• CLI Quickstart: run the bundled example and inspect the enrichment TSV.
• Expression20 vs GO: enrich the small expression-derived panel against official GOA human/GO targets.
• GO Impact Report: compare 2021 vs 2026 GO/GOA over 500 expression-derived gene sets using a parameterized report config, then inspect threshold crossings and timings.

CLI Quickstart

Source notebook: ../../notebooks/01_cli_quickstart.ipynb

Start with the common bioinformatics case: a one-column list of human genes from an over-expression result, and return a table of enriched GO terms. This example uses 50 interferon/inflammatory genes and a prepared official GOA human/GO snapshot.

Output:

Output:

mkdir -p notebooks/generated
genesets-rs run examples/human-go-enrich.yaml > notebooks/generated/quickstart_go.tsv
echo "Wrote notebooks/generated/quickstart_go.tsv"

Output:

Wrote notebooks/generated/quickstart_go.tsv

Output:

Output:

Broad GO terms are often expected at the very top. A useful first pass is to also inspect terms below a target-size threshold.

Output:

Expression20 vs GO

Source notebook: ../../notebooks/02_expression20_eval.ipynb

This notebook runs 20 expression-derived MyGeneset/MSigDB signatures against an official GOA human/GO snapshot. The CLI does the enrichment; Python only summarizes the Parquet results.

python3 scripts/fetch_mygeneset_eval.py \
  --manifest evals/expression20/sets.tsv \
  --out-dir evals/expression20/generated \
  --skip-existing

Output:

Using existing MyGeneset snapshot at evals/expression20/generated

python3 scripts/run_expression20_go_demo.py --skip-existing

Output:

reused_existing: notebooks/generated/expression20_go_now.parquet
metadata: notebooks/generated/expression20_go_now.yaml

Output:

Output:

Output:

The strongest unfiltered GO hits are often broad terms. For a more useful first look, the next table keeps terms with target size at most 500 genes.

Output:

Five-Year GO Impact Report

Source notebook: ../../notebooks/03_go_diff_with_duckdb.ipynb

This tutorial treats GO temporal change analysis as a configured convenience report. The report command runs the two enrichments, performs the threshold-crossing comparison, writes Parquet plus YAML/JSON metadata, and records the timing and parameters used. The notebook then inspects those report artifacts with DuckDB instead of rebuilding the workflow step by step.

uv run --project python/genesets-workflows genesets-workflows go-impact evals/go_impact_5y_expression500.yaml

Output:

Wrote notebooks/generated/go_impact_5y_expression500/summary.yaml
Diff rows: gained_significant=23378, lost_significant=954, shared_significant=2760
Timings: old=3.925s, new=5.114s, compare=0.174s, total=9.744s

Output:

Output:

Output:

Output:

Output:

For each high-impact expression signature, the term columns identify the GO term with the largest signed change in -log10(Bonferroni adjusted p) after applying the configured ranked-term scope. In this report, candidates are current GO terms in goslim_generic or descendants of a goslim_generic term, with at most 1,000 annotated genes in the relevant snapshot. Positive values mean stronger enrichment in 2026; negative values mean weaker enrichment. For gained/lost calls, the non-significant side is anchored at the report cutoff (0.05) because this report stores significant-only result tables.

This is an inclusion slim plus a size guard. The complementary antislim is the same idea with the go_term_filter join inverted, which is useful when the goal is to explicitly remove a known broad region of the ontology from a report.

Output:

The next table asks which GO terms changed status across the largest number of expression signatures. It reuses the configured goslim_generic descendant scope and adds a stricter display guard of at most 500 genes so the table stays biologically interpretable.

Output:

Output:

Notebook Workflows

The CLI should remain the primary interface, but notebooks are useful for eval review and result exploration.

Recommended Shape

Use notebooks as consumers of CLI output and report artifacts, not as the primary source of truth. Scripts and Rust code remain responsible for reproducible eval generation.

A typical workflow:

genesets-rs run evals/expression20/config.yaml

Then inspect evals/expression20/generated/results.tsv in a notebook with pandas, Polars, DuckDB, or plotting libraries.

For larger temporal reports, prefer a parameterized report config:

genesets-workflows go-impact evals/go_impact_5y_expression500.yaml

The notebook should show the report parameters and generated artifacts, then query the Parquet outputs.

The committed demonstrator notebooks live in notebooks/ and are linked from the Tutorials chapter:

• 01_cli_quickstart.ipynb;
• 02_expression20_eval.ipynb;
• 03_go_diff_with_duckdb.ipynb.

Python Notebook

Python notebooks are the most portable option for collaborators:

• read TSV outputs with pandas or Polars;
• plot p-value distributions;
• compare ranks across ontology releases;
• join results to source metadata.

The demonstrator notebooks should use %%bash for CLI commands and Python for analysis only. They should not expose Rust APIs.

Repository Policy

Commit small, deterministic notebooks only when they add durable explanation. Keep outputs cleared or tiny, and write generated data under notebooks/generated/, ignored eval generated/ directories, or external artifacts.

Web Explorer

The explorer is a local browser for workflow report bundles. It is not deployed yet, and it does not run enrichment itself. It reads summary.yaml, query metadata, GMT query genes, and Parquet result/diff files from one or more existing report directories.

Install or run with the optional explorer dependencies:

just browser

This opens the default set of report bundles, currently:

• the 5,000-set current GOA all-vs-IBA comparison;
• the 5,000-set 2021-vs-2026 GO/GOA comparison;
• the 4,313-set current GOA all-vs-no-contributes_to comparison, excluding GO-derived query sets.

To open one directly:

just browser-iba
just browser-go5y
just browser-contributes

The equivalent explicit command for the IBA bundle is:

uv run --project python/genesets-workflows --extra explorer \
  genesets-workflows explore notebooks/generated/go_iba_impact_expression5000_diverse

To open the default browser automatically:

uv run --project python/genesets-workflows --extra explorer \
  genesets-workflows explore notebooks/generated/go_iba_impact_expression5000_diverse --open

The server defaults to http://127.0.0.1:8765.

Bundle Inputs

A bundle is any directory with a summary.yaml written by genesets-workflows go-impact. The summary points to:

• queries.gmt;
• queries.metadata.json;
• A-target and B-target result Parquet files;
• threshold-crossing diff Parquet;
• optional term-coverage Parquet files.

If no bundle path is supplied, the explorer looks for notebooks/generated/*/summary.yaml.

For the artifact contract, see Report Bundles.

Current Views

The first screen ranks gene sets. The table can be searched and filtered by source family or stratum, then sorted by criteria such as:

• specific IBA losses;
• all IBA losses;
• specific IBA gains;
• largest p-value delta;
• result rows in either run;
• gene count.

Selecting a gene set opens detail panels:

• Diffs: lost, gained, and shared enrichment calls for the configured comparison;
• A targets: retained enriched targets from the left run;
• B targets: retained enriched targets from the right run;
• Genes: the query genes from the selected gene set.

The current gene panel shows query genes. To show per-target overlap gene lists, the workflow layer should add a target-membership artifact such as target_gene_membership.parquet and let the explorer join query genes to target genes on demand.

Design Boundary

The web explorer belongs to the workflow layer. It is a reader over report artifacts and DuckDB queries; it should not add web-specific behavior to the Rust enrichment kernel.

Relationship To Curated Gene Set Browsing

The current explorer browses eval report bundles. A curated gene set browser is a related but separate surface: it should read validated interpretation YAML or generated JSON from curation/, then render one page per curated gene set plus indexes by context, collection, category, and recovery status.

Both browsers can share workflow-layer tooling and static assets, but they should keep different input contracts:

• eval explorer: summary.yaml plus Parquet report artifacts;
• curated browser: LinkML-validated YAML plus generated corpus indexes.

See Static Browser Plan.

Report Bundles

A report bundle is a directory of immutable artifacts produced by a workflow run. Bundles are the contract between batch computation, notebook analysis, docs summaries, and the local web explorer.

The local explorer currently expects summary.yaml from genesets-workflows go-impact, plus the files referenced by that summary.

Minimal Shape

run-dir/
  summary.yaml
  summary.json
  queries.gmt
  queries.metadata.json
  left-results.parquet
  right-results.parquet
  left-vs-right.diff.parquet
  left-vs-right.diff.yaml

The exact filenames may vary by report. summary.yaml is the stable entry point and should record relative paths, source URLs, file digests, parameters, row counts, timing, and version labels.

Why Bundles

Bundles keep the interfaces simple:

• Rust writes result tables and comparison outputs.
• Python records source and report metadata.
• DuckDB queries Parquet without loading all rows into memory.
• Notebooks and docs can reproduce summaries from the same artifacts.
• The web explorer can browse a report without rerunning enrichment.

Generated bundles should usually live under an ignored directory such as notebooks/generated/ or an explicit external results directory. Small eval configs and manifests belong in git; large generated report artifacts usually do not.

Browser Requirements

For the current explorer, a bundle should provide:

• query genes as GMT;
• query metadata as JSON;
• left and right result Parquet files;
• threshold-crossing diff Parquet;
• optional term-coverage Parquet files.

Future bundle formats should preserve the same design: one small manifest file points to typed data artifacts, and consumers read the manifest instead of guessing file names.

Curated Gene Set Interpretations

The curation corpus stores expert GO interpretations of non-GO gene sets. It is both a precision/recall fixture for enrichment reports and the source content for a future curated gene set browser.

Interpretations live in curation/genesets/*.yaml. Each file conforms to the LinkML schema at curation/schema/genesets_interpretation.yaml and is validated with:

• LinkML structural validation;
• linkml-term-validator for ontology IDs and labels;
• linkml-reference-validator for cited literature snippets.

What A Curated File Contains

Each curated interpretation records:

• stable source identifier and source collection;
• biological context terms such as cell type, disease, phenotype, perturbagen, or experimental condition;
• curated GO term associations;
• curator role category, confidence, and specificity;
• recovery status, distinguishing annotation gaps from gene set membership gaps;
• enrichment stats that seeded an association;
• optional cited evidence and curator notes.

The current corpus includes MSigDB C8 cell-type signatures, KEGG legacy disease and pathway sets, Reactome and Hallmark sets, and a C7 perturbation contrast. The manifest at curation/genesets/manifest.tsv is the index of curated sets.

Why This Is Separate From Evals

Eval result tables say which GO terms a run recovered. The curation corpus says which GO terms are biologically appropriate for a gene set. Keeping those separate lets reports measure precision and recall without rewriting biological truth to match the current state of GO annotations.

category is the biological judgment. recovery_status explains whether the current annotation and membership state makes that judgment recoverable:

• annotation_supported: current annotations recover the term;
• annotation_gap: relevant genes are present, but GO annotation is too shallow;
• membership_gap: the gene set lacks genes needed to support the term.

Local Commands

Validate the corpus and build the current curation report:

just curate-validate
just curate-report

See Curation Workflow for the full workflow and Schema Reference for the LinkML model.

Browser Direction

The curated corpus should become browsable as generated pages or a static site. The source of truth should remain the validated YAML files, with generated JSON and HTML/Markdown as read models.

Use LinkML-generated docs for schema reference. Use generated per-gene-set pages for the corpus itself. See Static Browser Plan.

Curation Workflow

The curation workflow turns enrichment output into LinkML-validated biological interpretations of non-GO gene sets.

Local Setup

Use the curation extra from the Python workflow package:

uv run --project python/genesets-workflows --extra curation \
  genesets-workflows curate --help

The justfile provides the common repository commands:

just curate-validate-schema
just curate-validate
just curate-report

Steps

1. Validate the schema’s ontology-backed enum meanings:

   just curate-validate-schema
   

2. Draft a curated interpretation from enrichment output:

   uv run --project python/genesets-workflows --extra curation \
     genesets-workflows curate draft MSIGDB:<SET> \
     --enrichment-tsv path/to/enrichment.tsv \
     -o curation/genesets/<SET>.yaml
   

3. Adjudicate the draft by assigning category, confidence, specificity, and recovery_status. Add curator-only core terms and literature evidence where needed.

4. Validate every curated YAML file:

   just curate-validate
   

5. Build the precision/recall report:

   just curate-report
   

Judgment Model

category is the biological judgment. A term that is central to the gene set’s biology should remain core even if current GO annotations do not recover it.

recovery_status records whether the current annotations and membership make that biological judgment recoverable:

• annotation_supported: current annotations recover the term.
• annotation_gap: relevant genes are present, but GO annotation is too shallow.
• membership_gap: the gene set lacks genes needed to support the term.

This separation keeps the corpus useful both as an enrichment eval fixture and as a source of curation targets.

Schema Reference

The curation schema lives at curation/schema/genesets_interpretation.yaml. It is a LinkML schema for a single curated interpretation of a non-GO gene set.

Main Classes

• GeneSetInterpretation: one curated gene set interpretation.
• BiologicalContext: ontology terms describing what the set is about, such as cell type, disease, phenotype, perturbagen, or experimental condition.
• TermAssociation: one curated association between the gene set and a GO term.
• Term: an ontology term reference with an ID and label.
• EnrichmentStats: enrichment values that seeded an association.
• EvidenceItem: literature evidence supporting, refuting, or contextualizing a term association.

Important Fields

• gene_set_id: stable source identifier, such as MSIGDB:<SET>.
• contexts: ontology-grounded biological context for the gene set.
• associations.term: the GO term being judged.
• associations.category: curator role judgment.
• associations.confidence: confidence in the biological judgment.
• associations.specificity: how specific the term is to the context.
• associations.recovery_status: whether current annotations and membership make the term recoverable.
• associations.evidence: cited evidence with optional validated snippets.

Validation

Validation has three layers:

• LinkML structural validation checks required fields and enum values.
• linkml-term-validator checks ontology IDs and labels for GO, CL, UBERON, MONDO, CHEBI, PR, HP, NCBITaxon, EFO, and related prefixes.
• linkml-reference-validator checks cited evidence snippets against referenced literature when snippets are provided.

Use:

just curate-validate-schema
just curate-validate

LinkML Browser

A generated LinkML schema browser is useful as reference documentation for the schema. It should complement, not replace, a curated gene set browser. Most users will want pages organized by gene set and biological context; schema docs are organized by classes and slots.

Static Browser Plan

The curated interpretation corpus should become browsable from docs or a dedicated static site before it needs a deployed dynamic service.

Recommended Direction

Generate static pages from curation/genesets/*.yaml:

curation/genesets/*.yaml
  -> validate with LinkML, term validator, and reference validator
  -> materialize a small JSON index
  -> generate one page per curated gene set
  -> generate collection, context, category, and recovery-status index pages

The generator should live in the Python workflow package, for example as:

genesets-workflows curate build-pages \
  --dir curation/genesets \
  --out-dir docs/src/generated/curation

Those generated pages can be committed while the corpus is small, or built in CI for a static docs deploy once the corpus grows.

Page Shape

Each gene set page should show:

• source ID, name, collection, taxon, direction, and gene count;
• biological contexts with ontology IDs and labels;
• curated GO associations grouped by role category;
• confidence, specificity, and recovery status;
• enrichment stats and overlap genes when present;
• evidence references and curator notes;
• links to neighboring index pages by collection, context, GO term, category, and recovery status.

LinkML Browser Role

Use generated LinkML docs for the schema reference. Use generated corpus pages for the curated instances. A schema browser explains what TermAssociation means; an instance browser explains what MSIGDB:HAY_BONE_MARROW_DENDRITIC_CELL means.

Web UI Relationship

The current web explorer browses eval report bundles. A future curated gene set browser can reuse the same principles but should read curated corpus artifacts instead of report bundles:

• source of truth remains the validated YAML corpus;
• generated JSON gives the browser a stable read model;
• static pages support simple hosting and durable links;
• a richer client-side UI can be added later over the same JSON.

Do not add curated-browser behavior to the Rust compute engine. If enrichment is needed to refresh stats, run that through workflows and write artifacts back to the curation/report layer.

Reference

These pages provide background and design context that cuts across the compute, workflow, explorer, and curation surfaces.

• Ontology Prep: preparing ontology tables and closure.
• Competitive Landscape: related tools and design comparisons.
• Roadmap: planned compute, workflow, explorer, and curation work.

Ontology Prep

The core engine assumes ontology prep has already happened. That keeps enrichment fast and makes the closure policy explicit.

Recommended First Path: OAK

OAK, the Ontology Access Kit, is a strong first prep layer because it already handles many biomedical ontology concerns:

• local files and ontology selectors;
• CURIEs, labels, aliases, mappings, and obsolete terms;
• command-line and Python access;
• graph traversal over configured predicates;
• cached downloads for common ontology backends.

OAK’s CLI follows the pattern:

runoak --input HANDLE COMMAND [ARGS]

For biological ontologies, OAK supports predicate shortcuts for is_a and part_of:

runoak -i sqlite:obo:go ancestors GO:0006915 -p i,p

For labels, fill-table can populate a table of IDs with labels:

runoak -i sqlite:obo:go fill-table terms.tsv

For this project, the prep target is two TSV files:

terms.tsv:   term_id, name
closure.tsv: child, ancestor

The closure should be reflexive and should document which predicates were included.

When Horned-OWL Makes Sense

Horned-OWL is a Rust OWL library focused on performance and large ontologies. It can parse and manipulate OWL ontologies directly in Rust, and its published docs describe 20x to 40x speedups versus OWL API-based workflows for some large-ontology validation tasks.

It is attractive if we want a fully native prep path, especially for:

• parsing OWL/XML, RDF/XML, Functional Syntax, or Manchester Syntax;
• extracting SubClassOf and selected object-property restrictions;
• building relation closure without Python;
• shipping a single Rust binary for prep plus enrichment.

The tradeoff is scope. Correct OWL reasoning is more than graph walking. For MVP prep, OAK or precomputed relation-graph exports are lower risk. A future genesets-prepare crate can use Horned-OWL for the fast graph-walk case and clearly label it as such.

Closure Policy

Every prepared closure should record:

• ontology source and version;
• relation policy, such as is_a only or is_a + part_of;
• whether inferred axioms were included;
• whether obsolete terms were dropped, migrated, or retained;
• whether the closure is reflexive.

Those details affect enrichment results and must be part of eval metadata.

GAF Filtering Helpers

GAF filtering belongs in prep, not in the enrichment core. For GO evals, the helper script scripts/prepare_go_eval.py implements useful comparison variants:

• always remove NOT annotations;
• optionally remove contributes_to;
• restrict to IBA;
• restrict to IBA plus IEA.

The script reads GAF qualifier and evidence-code columns directly, writes ordinary gene_id, term_id annotation tables, and records the exact filter policy in YAML metadata.

Competitive Landscape

The goal is not to be the broadest enrichment product. It is to be a very fast, reproducible core engine with an eval framework.

Existing Tools And What They Optimize

GOATOOLS is a Python library and CLI for GO analysis. It parses OBO, GAF, GPAD, NCBI gene2go, and related formats; it propagates counts to parent terms by default; and it uses Fisher exact tests with many correction options.

g:Profiler and gprofiler2 provide a maintained web/API/R workflow across GO, pathways, disease annotations, and identifier namespaces. They use hypergeometric enrichment and offer g:SCS, Bonferroni, and FDR correction choices.

clusterProfiler is the dominant R/Bioconductor ecosystem tool for ORA and ranked-list workflows. It has GO-specific helpers and a universal enricher interface using user-provided TERM2GENE mappings.

topGO is important because it exposes topology-aware GO methods such as classic, elim, weight, weight01, lea, and parentchild with multiple statistics.

Ontologizer is important for parent-child and model-set ideas. Its papers frame the issue that GO terms are dependent through inheritance, so standard term-for-term Fisher tests can be biologically redundant.

Enrichr is a popular web/API ecosystem with many gene set libraries, fast interactive workflows, and p-value, z-score, combined-score, and adjusted p-value outputs.

PANTHER powers the GO Consortium enrichment service and is a strong reference point for user expectations around GO term enrichment and background selection.

MyGeneset.info is not a direct competitor. It is a useful public gene-set source and API for eval inputs because it exposes GO, MSigDB, Disease Ontology/OMIM, WikiPathways, Reactome, KEGG, and species filters.

Gene-Set-Informed Factor Models

A separate paradigm consumes gene sets without doing over-representation analysis at all. Instead of asking “is this set enriched in my ranked or thresholded gene list?”, these methods treat curated sets as priors on the latent structure of an expression matrix and let the model decide, per sample or per cell, which annotated programs are active.

f-scLVM (Buettner et al., 2017), implemented in the slalom package for R and Python, is the reference example. It is a sparse Bayesian factor analysis model for single-cell RNA-seq in which each annotated factor is initialized from a gene set (typically MSigDB or Reactome pathways): the set membership seeds the factor’s loading mask. The model then, jointly across all factors:

• estimates a per-factor relevance so it can switch annotated factors on or off rather than forcing every pathway to explain variance;
• refines the gene-set annotation itself, adding genes that load on a factor but were not in the input set and down-weighting members that do not;
• infers additional unannotated sparse and dense factors to absorb residual structure (technical and biological) that no input set explains.

The output is interpretable per-cell factor activations attached to named pathways, which is why it is used to identify and annotate cell subpopulations rather than to rank pathways for a single contrast.

expiMap (Lotfollahi et al., 2023), implemented in the scArches ecosystem, is the deep-learning successor to the same idea. It is a variational autoencoder whose latent nodes are gene programs: a binary mask ties the decoder so that each latent dimension reconstructs only the genes of its annotated set, making the otherwise opaque latent space directly interpretable. It carries over every move f-scLVM makes — it learns a per-cell activity for each program, soft-prunes uninformative programs (group-lasso regularization on the decoder, the relevance analog), and refines membership by allowing the mask to admit a few extra genes per program. Its headline addition is a set of add-on de novo learnable nodes: unconstrained latent dimensions that capture biology the annotations miss, the neural-network counterpart of f-scLVM’s unannotated factors. Because it is built on scArches, its distinctive use case is reference mapping — projecting a new query dataset onto an existing atlas’s gene-program latent space and reading off which programs differ, including in perturbation and disease-response settings.

The trade-off between the two mirrors linear-vs-nonlinear generally: f-scLVM is a linear factor model that is fast, fully Bayesian, and easy to reason about; expiMap is a nonlinear VAE that scales to atlas-sized data and supports transfer learning, at the cost of interpretability being mediated by the architecture rather than read directly off loadings.

This matters to a gene-set engine for two reasons. First, it is a downstream consumer of exactly the libraries this project scores and curates, so set quality propagates into factor interpretability. Second, the “refine the annotation” step is the same concern the curated corpus encodes by hand: the nominal membership of a named set and the genes that actually behave like the set are not identical. Both f-scLVM and expiMap learn that gap from data (membership refinement); the curation gold standard records it from expert judgment (recovery_status / membership_gap). The two are complementary views of annotation noise, not competitors.

GSEA and ssGSEA (single-sample GSEA) sit between ORA and factor models: they are rank-based set-scoring methods that, in the ssGSEA case, also produce a per-sample per-set activity score. They share the “score named programs per sample” goal with f-scLVM but do not jointly model factors, switch sets off, or refine membership.

Connection To The Curated Corpus

The curated gold standard in curation/ is, structurally, a hand-built version of what these models estimate. A single non-GO gene set is not curated to one label; it is decomposed into multiple GO programs, each with a curator judgment. HALLMARK_INTERFERON_GAMMA_RESPONSE, for example, resolves into a response-to- type-II-interferon program, an antigen-processing/MHC-class-I program, an antiviral-defense program, and a nonspecific translation residual. That is a factor decomposition of the set, recorded with evidence rather than fit from an expression matrix, and the schema axes line up with the model concepts:

• multiple associations per set ≈ the several latent program nodes one set seeds;
• category (core_process … nonspecific / false_association) ≈ per-factor relevance / group-lasso pruning — the nonspecific translation term is the factor f-scLVM would drive toward zero relevance;
• recovery_status: membership_gap ≈ data-driven membership refinement: expiMap’s reported addition of nine B-cell markers to its predefined BCR program is exactly a membership_gap correction, learned from data instead of curated;
• series / contrasting poles ≈ a program’s activation contrasting across cell states.

Two consequences. First, the corpus is natural supervision and a benchmark for the membership-refinement step specifically: the models propose “gene X belongs in program P”; the gold standard holds independent, cited judgments of that same claim. Second, the alignment is only with membership_gap — the models refine set membership but never touch the ontology, so annotation_gap (GO too shallow) stays orthogonal, just as it does for the enrichment tools above.

The overlap is concrete, not aspirational: f-scLVM seeds factors from MSigDB Hallmark and Reactome and headlines a G2/M-checkpoint cell-cycle factor; expiMap trains on Reactome plus PanglaoDB marker sets and headlines interferon programs. The corpus already curates that family — HALLMARK_G2M_CHECKPOINT, HALLMARK_INTERFERON_GAMMA_RESPONSE, and the HAY_BONE_MARROW_* / DESCARTES_* cell-type marker sets — so the priors these models consume are the same objects this project scores. Note the scope line, though: this repo curates gene sets and their GO interpretations, not expression matrices. The in-scope way to bring a signature from one of these papers into the corpus is the existing LIT:DISEASE_ACTIVITY pattern (a derived gene list with the paper as identity, membership, and evidence), not hosting the underlying single-cell data.

Our Wedge

The core differentiators should be:

• ontology-neutral input tables rather than GO-specific parsing in the hot path;
• fast dense bitset scoring for repeated jobs;
• first-class N x N term/set matrix mode;
• explicit background semantics;
• reproducible prep metadata;
• benchmark and correctness evals as part of the repo.

What Not To Compete On Yet

We should not try to match all visualization, identifier conversion, web UI, or database coverage of mature ecosystems in the core crate. Those are wrapper-layer concerns.

The MVP should instead make it easy to prove:

• this input produced this table;
• this table used this closure policy;
• these p-values match the reference implementation;
• this workload scored faster than comparable tools.

Research Sources

• GOATOOLS project page: https://pypi.org/project/goatools/
• gprofiler2 vignette: https://rdrr.io/cran/gprofiler2/f/inst/doc/gprofiler2.Rmd
• topGO manual: https://bioconductor.org/packages/devel/bioc/vignettes/topGO/inst/doc/topGO_manual.html
• Ontologizer paper: https://academic.oup.com/bioinformatics/article/24/14/1650/182451
• f-scLVM / slalom paper: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1334-8
• expiMap paper: https://www.nature.com/articles/s41556-022-01072-x
• Enrichr paper: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-14-128
• GO Consortium enrichment page: https://geneontology.org/docs/go-enrichment-analysis/
• MyGeneset.info docs: https://docs.mygeneset.info/
• OAK CLI docs: https://incatools.github.io/ontology-access-kit/cli.html
• Horned-OWL docs: https://docs.rs/horned-owl/latest/horned_owl/

Roadmap

Short Term

• Add Benjamini-Hochberg FDR.
• Add top-k and p-value-threshold output modes.
• Add serialized prepared indexes.
• Reframe docs around compute, workflows/evals, local explorer, and curated interpretation corpus.
• Add eval adapters for GOATOOLS and clusterProfiler.
• Add file digests and ontology metadata joins to compare metadata.
• Keep eval runners on Parquet by default.
• Move remaining ad hoc eval/source scripts into genesets-workflows.
• Add workflow config schema validation and genesets-rs version checks.
• Add a genesets-workflows curate build-pages command or equivalent static corpus page generator.

Medium Term

• Add OAK-based ontology prep examples.
• Add native Horned-OWL graph-walk prep as an optional crate or feature.
• Add ranked-list input schema.
• Add ranked-list enrichment statistics.
• Add richer benchmark reports for GitHub Pages.
• Add DuckDB notebooks over Parquet result and diff tables.
• Add hierarchy-aware Reactome prep and compare it with the flat baseline.
• Generate LinkML schema docs and curated gene set instance pages for the docs site.
• Add a curated gene set browser over generated corpus JSON.

Longer Term

• Add topology-aware model families.
• Explore gene-set-informed factor models (f-scLVM / slalom and expiMap-style) as downstream consumers of curated libraries, and reuse the curated recovery_status / membership_gap judgments as a benchmark for data-driven annotation refinement. See the Competitive Landscape.
• Add richer version-diff reports for ontology release comparisons.
• Add wrapper layers for species, release metadata, and identifier mapping.
• Add visualization notebooks that consume normalized result TSVs.
• Deploy the eval explorer and curated corpus browser from stable generated artifacts.