sharp

CLAUDE.md — S(H)ARP Project Context

Read this before touching any code. Additional detail in docs/.

What this project is

S(H)ARP predicts Biosynthetic Gene Clusters (BGCs) in Streptomyces and related actinomycetes using SARP transcription factors as anchors. A SARP (Streptomyces Antibiotic Regulatory Protein) is a regulator with an HTH-BTAD domain (± NB-ARC / TPR / AAA / LuxR) that binds heptameric repeats (afsR-box) in BGC promoters.

Differentiator vs. antiSMASH: we use regulatory context + protein language model embeddings, not just biosynthetic enzyme patterns.

See docs/PIPELINE.md for the full biological pipeline. See docs/ARCHITECTURE.md for module ownership.

Repo layout

project_root/
├── CLAUDE.md                   ← you are here
├── .env.example                ← template for machine-specific settings (copy → .env, gitignored)
├── pyproject.toml              ← editable install: `pip install -e .`
├── pixi.toml                   ← environment (use pixi, not conda/pip directly)
├── src/sharp/
│   ├── __init__.py
│   ├── config.py               ← paths + `.env` loading + config dataclasses (DONE)
│   ├── io.py                   ← all data types + file I/O (DONE)
│   ├── metrics.py              ← pure metric math (DONE)
│   └── evaluate.py             ← benchmark step orchestration (DONE)
├── scripts/
│   ├── generate_mock_data.py           ← synthetic proteins for embedding step (DONE)
│   ├── generate_mock_benchmark_data.py ← synthetic predictions + GT for benchmark (DONE)
│   └── prepare_mibig_ground_truth.py   ← parse MiBIG 4.0 JSON → ground_truth.tsv (DONE)
├── tests/
│   ├── conftest.py
│   ├── test_config.py
│   ├── test_io.py
│   ├── test_model_management.py
│   ├── test_extract_embeddings.py
│   ├── test_generate_mock_data.py
│   ├── test_metrics.py
│   ├── test_evaluate.py
│   └── test_prepare_mibig.py
└── data/
    ├── raw/          ← immutable inputs (MiBIG dump, downloaded genomes)
    ├── interim/      ← intermediate pipeline artifacts
    ├── processed/    ← final outputs (model.pkl, report.html, benchmark.json)
    └── mock/         ← synthetic data for testing

Note: extract_embeddings.py and model_management.py were implemented and tested but belong in src/sharp/ — they may already be there if you’ve been working in this repo. If missing, see docs/ARCHITECTURE.md for their specs.

Git workflow

Commit after every meaningful unit of work. Follow Conventional Commits:

<type>(<scope>): <short description>

[optional body]

Types used in this project:

Type When
feat new pipeline step, new script, new metric
fix bug fix in existing code
test adding or fixing tests
refactor restructuring without behavior change
docs CLAUDE.md, docs/, docstrings
chore pixi.toml, pyproject.toml, CI config
data scripts that produce or transform data files

Scopes are module names or script names: io, metrics, evaluate, extract-embeddings, prepare-mibig, run-antismash, etc.

Examples:

feat(metrics): add reciprocal_overlap and BenchmarkResult
test(evaluate): add end-to-end orchestration tests
fix(io): handle missing region_id in FASTA header gracefully
feat(prepare-mibig): add --inspect mode for schema verification
data(prepare-mibig): build Streptomyces ground truth from MiBIG 4.0
refactor(extract-embeddings): extract residue_mean_pool as pure function
docs(claude): add benchmark comparison section and backlog tier 0
chore: add antismash and deepbgc to pixi.toml

Keep commits atomic — one logical change per commit. Don’t batch unrelated changes (e.g. don’t fix a bug and add a feature in the same commit).

Environment

pixi run python ...          # always use pixi, not bare python
pixi run pytest              # run tests
pixi run python -m sharp.evaluate --help

Package is installed editable: import sharp.io works anywhere.

Conventions — read before writing any code

Coordinates: 0-based half-open [start, end) everywhere. MiBIG uses 1-based inclusive [from, to]; prepare_mibig_ground_truth.py converts on ingest (start - 1, end unchanged). Never store 1-based coords in any data type or file.

Data types live in io.py, not in a separate types.py. The rule: a type lives in the module that constructs it. Extract to types.py only if a third module needs it without going through io.

Config dataclasses in config.py. One frozen dataclass per pipeline step (e.g. EmbeddingConfig, EvaluateConfig). Steps receive a config object, not loose **kwargs.

.env is for machine identity, not pipeline behaviour. config.py reads an optional gitignored .env at the project root (template: .env.example) to resolve the data-dir constants and DEFAULT_DEVICE. Precedence: real env var → .env → in-repo default.

Belongs in .env Never in .env
SHARP_DATA_ROOT (+ SHARP_{RAW,INTERIM,PROCESSED,MOCK}_DIR) Thresholds (min_cluster_frac, min_p_bgc, e-value cutoffs)
SHARP_DEVICE — hardware physically present model_name, batch_size, max_length
Credentials, once a step actually needs one Anything that changes what a run means

Reason: a benchmark number must be reproducible from its command line. If a threshold could come from an unversioned file, the same command could produce different numbers on the laptop and the server. When adding a new step, put its knobs in the config dataclass and on the CLI — reach for .env only if the value describes the machine rather than the experiment.

Paths are resolved once at import, so tests that need different values reload the module under a patched environment (see tests/test_config.py). Don’t add python-dotenv — the loader in config.py is intentionally ~15 lines.

Each pipeline step = one module with a run(cfg: StepConfig) -> None function and a build_parser() -> argparse.ArgumentParser function. Entry point: python -m sharp.<step>.

When you add a new file (script or module), do these two things in the same change:

  1. Update the directory-structure tree in README.md (## Directory Structure) so it stays accurate.
  2. Document its usage in both README.md (a runnable command block, like the “Preparing MiBiG / BGC Atlas Database” sections) and CLAUDE.md. The single source of truth is the file’s own module docstring — mirror its Usage: block; don’t invent new invocations. Keep the three (docstring, README, CLAUDE.md) in sync.

Side-effect isolation: metrics.py is pure (no I/O, no logging). io.py owns disk. Orchestration modules (evaluate.py, extract_embeddings.py, etc.) call both and log.

Streaming writes for large files. Parquet is written batch-by-batch via pq.ParquetWriter context manager. Never accumulate all rows in memory.

Tests mirror src/sharp/ one-to-one. Test file for sharp/foo.py → tests/test_foo.py. Scripts tested in tests/test_<script_name>.py with sys.path injection (see existing examples).

Monkeypatching rule: patch on the importing module, not the source. If extract_embeddings.py does from sharp.model_management import Embedder, patch sharp.extract_embeddings.Embedder, not sharp.model_management.Embedder.

What is DONE (with tests)

Module / Script Responsibility Tests
sharp/config.py Paths (env-overridable), .env loading, DEFAULT_DEVICE, EmbeddingConfig, EvaluateConfig test_config.py
sharp/io.py ProteinRecord, PredictedRegion, KnownCluster; FASTA r/w; parquet r/w; TSV r/w; JSON w test_io.py
sharp/model_management.py ESM-2 registry, device selection, residue_mean_pool, Embedder, ensure_model_available test_model_management.py
sharp/extract_embeddings.py Embedding extraction step: load FASTA → embed → write parquet test_extract_embeddings.py
sharp/metrics.py overlap_bp, merge_intervals, covered_bp, matches, reciprocal_overlap, MatchCriterion, evaluate_predictions, BenchmarkResult test_metrics.py
sharp/evaluate.py Benchmark step: load predictions + GT (+ optional --contigs scope) → compute metrics → write JSON test_evaluate.py
scripts/generate_mock_data.py Synthetic proteins → FASTA (for embedding step smoke tests) test_generate_mock_data.py
scripts/generate_mock_benchmark_data.py Synthetic predictions + GT with controlled overlap (for benchmark smoke tests) test_evaluate.py (integration)
scripts/prepare_mibig_ground_truth.py MiBIG 4.0 JSON dir → ground_truth.tsv; handles 3.x fallback; --inspect mode; --genus and --exclude-eukaryotes taxonomic scoping; rejects unusable accessions (protein WP_/NP_, assembly GCA_/ASMnnnvn, WGS master ...01000000) and spans <500 bp, and collapses duplicate loci filed under two cluster ids test_prepare_mibig.py
scripts/prepare_bgcatlas_ground_truth.py BGC Atlas .gbk dump → bgcatlas_ground_truth.tsv (secondary/noisy GT) test_prepare_bgcatlas.py
scripts/convert_antismash_to_parquet.py antiSMASH sequence.json → predictions.parquet; no coord conversion; --inspect mode test_convert_antismash.py
scripts/convert_deepbgc_to_parquet.py DeepBGC .bgc.tsv → predictions.parquet; no coord conversion; --inspect mode test_convert_deepbgc.py
scripts/convert_gecco_to_parquet.py GECCO .clusters.tsv → predictions.parquet; start-1 coord conversion; --inspect mode test_convert_gecco.py
scripts/run_antismash.sbatch Slurm job: antiSMASH on the benchmark genome. CPU-only like DeepBGC. It takes --cpus and hands it to its own module scheduler, but that scheduler parallelises very little, so a wide allocation is wasted. Sizing measured (seff on array job 45315 index 1, which ran this script’s own default genome AL645882.2: 5.40% CPU efficiency of 16 cores, 1.62 GB peak) → 4 cores / 4G, matching run_antismash_array.sbatch; walltime stays 12h since $1 may be an arbitrary genome. Paths are explicit at the top of the file shell, no test
scripts/run_deepbgc.sbatch Slurm job: DeepBGC on the benchmark genome. CPU-only on a GPU-free node (the python=3.7 env predates CUDA-capable TF). Sized from seff on job 42995: 0.86 cores, 1.71 GB, 32 min → 2 cores / 8G / 2h. hmmscan dominates the runtime but DeepBGC does not thread it, so the pipeline is serial. Paths are explicit at the top of the file shell, no test
scripts/run_antismash_array.sbatch Slurm job array: antiSMASH over the whole benchmark set, one task per genome (index → line of analyzed_contigs.txt). Resumable (skips genomes with existing output), per-genome output dir. Sizing measured (seff on job 45315, indices 1–2, 2026-08-19): ~3 min wall, 5.4%/7.8% CPU efficiency of 16 cores (~1.3 cores of real parallelism), 1.6 GB peak → 4 cores / 4G (walltime 4h, deliberate slack over the measured ~3 min). antiSMASH ignores most of --cpus, the same way DeepBGC did (sized 8, measured 0.86). Takes the genome list as $1, an index offset as $2 (line = SLURM_ARRAY_TASK_ID + OFFSET, offset defaults to 0) the genome directory as $3 (default data/raw/genomes, so a genome database staged by build_genome_manifest.py --link-dir reuses the same script with an assembly-keyed pool) and the output pool root as $4 (default <tool>/out_benchmark; a separate campaign needs its own pool because merge_predictions.py converts every directory under --input-dir): Slurm’s MaxArraySize caps the highest legal array index at 1000, so a pool over that is submitted as successive windows reusing indices 1..1000 over the one list, rather than slicing it into per-chunk files. The list is read per task at task start, so editing it mid-flight renumbers lines under pending tasks. Logs go to logs/ (Slurm will not create it). shell, no test
scripts/run_deepbgc_array.sbatch Slurm job array: DeepBGC over the benchmark set, same shape. Per-task sizing carried from the measured single-genome run (2 cores / 8G; walltime 6h, deliberate slack over the measured 32 min); DeepBGC is single-threaded so the per-task shape does not change with the array. Takes the genome list as $1, an index offset as $2 (line = SLURM_ARRAY_TASK_ID + OFFSET, offset defaults to 0) the genome directory as $3 (default data/raw/genomes, so a genome database staged by build_genome_manifest.py --link-dir reuses the same script with an assembly-keyed pool) and the output pool root as $4 (default <tool>/out_benchmark; a separate campaign needs its own pool because merge_predictions.py converts every directory under --input-dir): Slurm’s MaxArraySize caps the highest legal array index at 1000, so a pool over that is submitted as successive windows reusing indices 1..1000 over the one list, rather than slicing it into per-chunk files. The list is read per task at task start, so editing it mid-flight renumbers lines under pending tasks. Logs go to logs/ (Slurm will not create it). shell, no test
scripts/run_benchmark.sh Step 4 of a benchmark run: merge the shared pool down to one scope, then evaluate. Derives all paths from a single scope name (data/interim/<scope>/{analyzed_contigs.txt,benchmark_ground_truth.tsv} → <tool>_predictions_<scope>.parquet → benchmark_<scope>_<tool>[.<label>].json), which makes the scope/GT mismatch structurally impossible. Pool at $POOL_ROOT/<tool>/out_benchmark (default ~/projects). Refuses to overwrite (--force), warns when the pool is smaller than the scope, reuses a parquet newer than the pool (--remerge forces) so --min-p-bgc sweeps are cheap; args after -- forward to sharp.evaluate. Multi-user — several people benchmark inside one server clone, so: --label NAME gives each person (or sweep point) its own JSON, the overwrite refusal names the file’s owner, the shared parquet is written to a temp path and renamed so a concurrent reader never sees a partial write, a provenance block (user, host, time, git commit, pool + genome count, evaluate args) is stamped into every result, and both output dirs are checked for group-writability up front. Does not select, download, or submit sbatch — those stay manual shell, no test
scripts/merge_predictions.py Per-genome baseline outputs → one predictions.parquet. Reuses each convert_<tool>_to_parquet.py (no reparsing — format assumptions stay in the tested converters); --tool antismash\|deepbgc\|gecco; --contigs reports genomes that produced no output at all, which would otherwise be scored as “tool found nothing”; --manifest genomes.tsv is required when the pool is keyed by assembly (a genome database) rather than by contig accession, or that check compares contigs against assembly directory names and flags every genome missing; --inspect mode test_merge_predictions.py
scripts/download_benchmark_genomes.sh Batch downloader: benchmark_genomes.tsv → data/raw/genomes/<ACC>.fasta, one per genome. Resumable (skips valid existing files, re-fetches truncated ones), retries with backoff, and verifies each FASTA header equals the expected accession — a mismatch would silently score that genome zero. Fetches by nucleotide accession, not from the cluster’s NCBI mirror (that mirror is assembly-indexed while the GT is nucleotide-keyed, and no index file bridges the two — see docs/NCBI_MIRROR.md) shell, no test
scripts/_fetch_nuccore.sh Sourced helper (not executed): fetch_nuccore (efetch + retry + “is this really FASTA” validation) and fasta_contig_ids. Shared by both downloaders so the rules live in one place shell, no test
scripts/download_genome.sh NCBI nuccore accession → data/raw/<ACC>.fasta + data/interim/analyzed_contigs.txt; defaults to AL645882.2 shell, no test
scripts/select_benchmark_genomes.py Ground truth → benchmark genome set: drops BGC-only deposits (--min-length), merges RefSeq/GenBank twins, ranks, and emits benchmark_genomes.tsv + analyzed_contigs.txt + benchmark_ground_truth.tsv (contigs normalized onto the primary accession). NCBI lengths cached in data/interim/record_lengths.tsv; --inspect/--offline modes test_select_benchmark_genomes.py
scripts/parquet_to_tsv.py Generic dump: any pipeline parquet file → TSV; list-typed columns (e.g. embeddings.parquet’s embedding vector) comma-joined per cell; --inspect mode test_parquet_to_tsv.py
scripts/build_genome_manifest.py Genome database (one NCBI assembly dump directory per genome) → genomes.tsv (assembly, contig, length, description), assemblies.tsv, assemblies.txt (job-array line list), analyzed_contigs.txt (scope file). Parses only the .fna (headers + streamed lengths; the .gbff is 2-5x larger and carries nothing needed). --link-dir stages <assembly>.fasta symlinks so the array scripts run over the database unchanged. Key is the GCF/GCA accession, not the directory name. --inspect/--limit modes test_build_genome_manifest.py
scripts/summarize_predictions.py Raw tool-vs-tool comparison over a genome database when there is no ground truth — no recall, no precision. Per (tool, threshold): n_regions, bp_called (merged intervals), regions_per_Mb, frac_bp_called, median_region_bp; plus a per-assembly table. Reports bp as well as counts because antiSMASH merges protoclusters while DeepBGC splits them, so counts are not commensurable; sweeps --thresholds because a scored tool’s count is a knob. Warns loudly when predictions do not join the manifest (dropped version suffix) test_summarize_predictions.py
scripts/predictions_to_ground_truth.py One tool’s predictions.parquet → ground-truth-shaped TSV, so evaluate.py measures tool-vs-tool agreement with no new metric code (coords pass through — both types are 0-based half-open). The resulting JSON’s recall means agreement with a tool, not with truth; run both directions test_predictions_to_ground_truth.py

What is NOT YET IMPLEMENTED

Implement these in order. Each is a pipeline step; each gets its own module + config dataclass + tests.

1. sharp/annotate.py — genome annotation

Input: data/raw/genome.fasta Output: data/interim/proteins.faa, data/interim/annotated.gbk, data/interim/genes.gff Tool: Bakta (shell out via subprocess) Config: AnnotateConfig(input_path, output_dir, threads, min_contig_length) Notes: Bakta writes its own output dir. Wrapper should copy/symlink the three output files to canonical interim paths. Validate that all three files exist after run.

2. sharp/detect_sarp.py — SARP detection by HMM

Input: data/interim/proteins.faa, data/raw/sarp_models.hmm Output: data/interim/anchors_sarp.tsv (columns: protein_id, contig, start, end, strand, score, type) Tool: hmmscan (shell out) Logic: Parse hmmscan tblout format. Add type=SARP. Filter by e-value threshold (default 1e-5). Config: DetectSarpConfig(proteins_path, hmm_path, output_path, evalue_threshold)

3. sharp/detect_heptarepeats.py — motif search in DNA

Input: data/raw/genome.fasta, FIMO motif file (afsR-box PWM) Output: data/interim/anchors_heptarepeats.tsv (same columns as above, type=heptarepeat) Tool: FIMO from MEME suite (shell out) Logic: Parse FIMO TSV output. Coords are already 0-based in FIMO output — verify this on real output before assuming. Config: DetectHeptarepeatsConfig(genome_path, motif_path, output_path, pvalue_threshold)

4. sharp/merge_anchors.py — unify anchor tables

Input: anchors_sarp.tsv, anchors_heptarepeats.tsv Output: data/interim/anchors.tsv Logic: Concatenate, deduplicate by position, sort by contig+start. Pure function, minimal I/O.

5. sharp/extract_neighborhood.py — genomic window extraction

Input: data/interim/anchors.tsv, data/interim/annotated.gbk Output: data/interim/neighborhoods.tsv, data/interim/neighborhood_proteins.faa, data/interim/neighborhood_dna.fna Logic:

6. sharp/annotate_domains.py — Pfam domain annotation

Input: data/interim/neighborhood_proteins.faa, data/raw/pfam_models.hmm Output: data/interim/domains.tsv (columns: protein_id, region_id, domain, e_value, start, end) Tool: hmmscan (shell out, domtblout format) Logic: Parse domtblout. One row per domain hit per protein. Filter by e-value. Note: region_id must be recovered from the FASTA header (use parse_fasta then build a protein_id → region_id map). Config: AnnotateDomainsConfig(proteins_path, hmm_path, output_path, evalue_threshold)

7. sharp/extract_kg_features.py — knowledge graph context features ⭐ yours

Input: data/interim/neighborhoods.tsv, data/interim/domains.tsv, data/raw/kg.gpickle Output: data/interim/kg_features.parquet (columns: region_id, n_similar_clusters, modal_class, has_large_sarp, ...) Logic: For each region, query the KG for clusters with similar domain architecture. Extract tabular features. Config: KgFeaturesConfig(neighborhoods_path, domains_path, kg_path, output_path) Note: KG is built by a separate one-time script (scripts/build_kg.py) — see docs/PIPELINE.md.

8. sharp/train.py — classifier training ⭐ yours

Input: data/interim/embeddings.parquet, data/interim/domains.tsv, data/interim/kg_features.parquet, data/raw/mibig_ground_truth.tsv Output: data/processed/model.pkl, data/processed/metrics.json, data/processed/feature_importance.tsv Logic:

9. sharp/predict.py — inference

Input: data/processed/model.pkl, same features as train Output: data/interim/predictions.parquet (columns: region_id, contig, start, end, p_bgc, predicted_class) Logic: Load model, run feature pipeline (same aggregation as train), predict. Output format must match PredictedRegion schema in io.py. Config: PredictConfig(...)

10. sharp/filter.py — heuristic post-filter

Input: data/interim/predictions.parquet, data/interim/domains.tsv Output: data/interim/filtered_predictions.parquet Logic (rules to start with):

11. sharp/generate_report.py — HTML report

Input: data/interim/filtered_predictions.parquet, data/interim/neighborhoods.tsv, data/interim/domains.tsv Output: data/processed/report.html Tool: Jinja2 Logic: One section per predicted BGC: region coordinates, class, p_bgc score, domain architecture diagram (SVG or simple HTML table).

Benchmark comparison — competitor baselines

Priority: high. The team wants S(H)ARP benchmarked against antiSMASH and DeepBGC (at minimum). Check recent literature for others.

The architecture is already correct: any tool’s output can be converted to predictions.parquet and passed through evaluate.py unchanged. Each tool gets one conversion script in scripts/.

Ground truth sources

Source Reliability Use as GT
MiBIG 4.0 ✅ Manually curated Primary — always use
BGC Atlas ⚠️ Computationally predicted, no manual curation Secondary — noisier, interpret separately

Benchmark scope caveat (verified 2026-07-29). Recall is measured only over ground-truth clusters on contigs the tool was actually run on — evaluate.py takes --contigs for this. Ground truth spans a database while a run spans one assembly, so without scoping, recall is capped by the ratio between them: on a real run of all three baselines against streptomyces_ground_truth.tsv (430 clusters / 363 contigs), the analyzed contig AL589148.1 carried exactly one coordinate-resolved cluster, and every tool scored recall ≤ 0.002. Pass the same --contigs file to every tool in a comparison — omitted, the scope is inferred from the predictions, which is optimistic (a contig analyzed but not called on drops out of the denominator) and logs a warning.

BGC Atlas results should be reported with a caveat in any paper/presentation: benchmark numbers on BGC Atlas are optimistic by nature (the positive labels are themselves predictions, so agreement with them doesn’t prove correctness).

MiBIG 4.0 coordinate-coverage caveat (verified 2026-07-07). ~45% of all MiBIG 4.0 entries — and 478 of 905 (53%) of Streptomyces entries — store their locus as location: {from: 0, to: 0}, i.e. the compound is characterized but the genomic coordinates are unknown. prepare_mibig_ground_truth.py correctly drops these (a coordinate-based benchmark can’t score a cluster with no interval; the drop count is logged as “N entries had no locus with usable coordinates”). The resulting Streptomyces ground truth is ~430 loci from 427 clusters, not ~900. Two consequences to report in any paper/presentation:

Baseline integration — converters, not wrappers

S(H)ARP never invokes the baseline tools. antiSMASH, DeepBGC, and GECCO each install into their own isolated pixi env under ~/.local/src/<tool>/ (via scripts/setup_<tool>.sh) — they have mutually incompatible dependencies and must stay isolated. You run each tool yourself (its own env, or HPC, or a container); S(H)ARP only parses the output files it leaves behind.

Setup scripts and .env. All three setup_<tool>.sh scripts source scripts/_load_env.sh (sourced, not executed) to read their write locations from .env: TOOLS_INSTALL_DIR (where the per-tool pixi envs go, formerly a hardcoded ~/.local/src repeated in all three) and, for antiSMASH/DeepBGC, ANTISMASH_DOWNLOADS_DIR / DEEPBGC_DOWNLOADS_DIR for the ~10GB / ~3GB reference databases. The helper uses set -a so the values are exported — DeepBGC has no path flag and reads DEEPBGC_DOWNLOADS_DIR from the environment itself; each script also applies a : "${VAR:=default}" fallback so it works with no .env. Note the README’s cd ~/.local/src/<tool> examples assume the default TOOLS_INSTALL_DIR; the scripts’ closing echos interpolate the real value. These two keys are the exception to the “no shell expansion” rule for .env: they use ${DATABASES:-$HOME/.local/share}, which bash resolves when sourcing but config.py’s minimal parser does not (os.path.expandvars has no :- support, so Python sees the literal string — inert, since no Python step reads them). DATABASES is intentionally not defined in .env.example: the server exports it, a laptop falls back to ~/.local/share, so neither machine needs an edit. Keep shell-only keys in that block and SHARP_* keys — read by Python — above it.

So each baseline gets one converter script (not a subprocess wrapper):

scripts/convert_<tool>_to_parquet.py --input <tool output> --output <predictions.parquet>

Runs entirely in the S(H)ARP env, no external binary, no tool-path config. Each converter isolates every tool-format assumption (column names, coordinate base) in one clearly-marked block and provides an --inspect mode that prints a real output file’s structure — verify the schema against actual output before trusting the parser (same pattern as prepare_mibig_ground_truth.py).

Coordinate base is tool-specific — verified per tool 2026-07-15 against a real run (antismash 8.0.4, deepbgc, gecco 0.10.3 on the same input FASTA, AL589148.1).

Tool p_bgc source Coordinate base Conversion
antiSMASH none → set 1.0 0-based half-open (verified) none
DeepBGC deepbgc_score 0-based half-open (verified) none — refutes old hypothesis
GECCO average_p 1-based inclusive (verified) start - 1, end unchanged — confirms old hypothesis

Evidence (span = end - start from the TSV/JSON row, cross-checked against the matching region/cluster .gbk LOCUS bp length, across every row in each output — not just one — since a single row can’t distinguish the two conventions if mis-signed. If span == LOCUS bp, the source is 0-based half-open. If span == LOCUS bp - 1, the source is 1-based inclusive (the .gbk extraction naturally has LOCUS bp = end - start + 1 bases for an inclusive interval):

Caution for implementation: a coordinate check on a single row can look consistent with either convention if you only compare one direction of the off-by-one; always check the full span-vs-LOCUS relationship (== vs == -1) across multiple rows before trusting the parser, as done above.

Tests parse a small, checked-in, real (trimmed) output fixture per tool — no tool execution in the suite (that would break env isolation). Same approach as test_prepare_mibig.py.

scripts/prepare_bgcatlas_ground_truth.py ✅ done (2026-07-07) Parses the BGC Atlas complete-bgcs dump — 204,661 antiSMASH-produced .gbk files, one region per file (downloaded by scripts/download_bgc-atlas.sh, DVC-managed under data/raw/complete-bgcs/). Output: data/raw/bgcatlas_ground_truth.tsv (same schema as mibig_ground_truth.tsv). Verified schema facts:

scripts/convert_antismash_to_parquet.py ✅ written (2026-07-15, verified against a real antismash 8.0.4 run; tests: tests/test_convert_antismash.py, fixture: tests/fixtures/antismash_sequence.json — trimmed real summary JSON) Parses sequence.json → data/interim/antismash_predictions.parquet. Iterate data['records'][*]['features'] where type == "region". Per region:

scripts/convert_deepbgc_to_parquet.py ✅ written (2026-07-15, verified against a real DeepBGC 0.1.0 run; tests: tests/test_convert_deepbgc.py, fixture: tests/fixtures/deepbgc_out.bgc.tsv — real, unmodified output) Parses out.bgc.tsv → data/interim/deepbgc_predictions.parquet. Real header (28 columns) confirms product_class exists as documented, but it sits at column 18, after several unrelated columns (detector_version, num_proteins, product_activity, per-activity probabilities) — don’t assume column order. Per row:

scripts/convert_gecco_to_parquet.py ✅ written (2026-07-15, verified against a real gecco 0.10.3 run; tests: tests/test_convert_gecco.py, fixture: tests/fixtures/gecco_sequence.clusters.tsv — real, unmodified output) Parses sequence.clusters.tsv → data/interim/gecco_predictions.parquet. Real header confirms sequence_id, cluster_id, start, end, average_p, max_p, type, plus per-class probability columns (nrp_probability, polyketide_probability, etc.) and proteins/domains list columns. Per row:

Running a full comparison

Run once, slice many. Baselines run over the broadest genome set once; every benchmark after that is a re-scope, not a re-run. The array output pool is keyed by accession (~/projects/<tool>/out_benchmark/<ACCESSION>/) and --contigs filters both the predictions and the ground-truth denominator (metrics.py), so a scope file plus its matching normalized GT fully define a benchmark. Never partition the pool per experiment — shared accession keys are what make re-slicing free. --contigs is therefore mandatory on merge_predictions.py, which would otherwise sweep in genomes from unrelated experiments. Name derived artifacts per scope (<tool>_predictions_<name>.parquet, benchmark_<name>_<tool>.json). See README → “Run once, slice many” and TODO.md for the queued scopes. docs/BENCHMARK_SCOPES.md is the catalogue: the pool (all 1,280 coordinate-resolved bacterial clusters / 1,112 accessions, BGC-only deposits included so they can be sliced out at scope time), every scope carved from it (_strep 156, _actino 575, _bact, per-genus negative controls, per-class slices, a 3-genome smoke scope), the rationale for each, and the exact commands that build the ground truths, the pool list, the download and the scopes.

The scaled run (50 genomes, 113 clusters) is the current published set — a single-genome run caps the recall denominator at 16 clusters. The single-genome flow below still works for a smoke test.

Scoreable MiBIG is much smaller than its entry count (measured 2026-08-21): 3,013 entries → 1,634 coordinate-resolved (1,420 accessions) → 1,280 bacterial (1,112) → 414 Streptomyces (352). After --min-length drops BGC-only deposits, full Streptomyces is 93 genomes / 156 clusters. --genus and --exclude-eukaryotes on prepare_mibig_ground_truth.py select these scopes.

# 1. Select the genome set from the ground truth. This is not just "sort by
#    cluster count": ~58% of MiBIG records are BGC-only deposits (the record IS
#    the cluster, so every tool scores ~1.0 by construction) and one physical
#    sequence can carry several accessions (NC_003888.3 and AL645882.2 are the
#    same S. coelicolor chromosome, with 15 clusters filed under one and 1 under
#    the other). Emits the set, the scope file, and a contig-normalized GT.
python scripts/select_benchmark_genomes.py \
    --ground-truth data/raw/streptomyces_ground_truth.tsv \
    --output-dir data/interim/benchmark_set

# 2. Fetch them (resumable, ~420 Mb). By nucleotide accession, so the FASTA
#    header IS the name the ground truth uses.
scripts/download_benchmark_genomes.sh

# 3. Run each baseline as a job array, then merge (see README), and evaluate
#    against benchmark_ground_truth.tsv + analyzed_contigs.txt from step 1 —
#    NOT the raw MiBIG ground truth, whose contig names are not normalized.

# ── single-genome smoke test ────────────────────────────────────────────────
# Fetch one genome and its --contigs scope file in one step.
scripts/download_genome.sh

# Build ground truth
python scripts/prepare_mibig_ground_truth.py \
    --input-dir data/raw/mibig_json_4.0 \
    --output data/raw/mibig_ground_truth.tsv --genus Streptomyces

# Run each baseline yourself in its own env (see scripts/setup_<tool>.sh), then
# convert its output — S(H)ARP never invokes the tools:
python scripts/convert_antismash_to_parquet.py \
    --input <antismash output dir/json> \
    --output data/interim/antismash_predictions.parquet

python scripts/convert_deepbgc_to_parquet.py \
    --input <deepbgc .bgc.tsv> \
    --output data/interim/deepbgc_predictions.parquet

python scripts/convert_gecco_to_parquet.py \
    --input <gecco .clusters.tsv> \
    --output data/interim/gecco_predictions.parquet

# Evaluate all against the same ground truth AND the same scope.
# --contigs lists the contigs the tools were run on (one per line, or a .fai);
# every tool must get the same file or the recall denominators differ.
# download_genome.sh already wrote analyzed_contigs.txt. For a genome obtained
# some other way, derive it the same way:
#   grep '^>' <genome.fasta> | cut -c2- | cut -d' ' -f1 > data/interim/analyzed_contigs.txt

for tool in antismash deepbgc gecco; do
    python -m sharp.evaluate \
        --predictions data/interim/${tool}_predictions.parquet \
        --ground-truth data/raw/mibig_ground_truth.tsv \
        --contigs data/interim/analyzed_contigs.txt \
        --output data/processed/benchmark_${tool}.json
done

python -m sharp.evaluate \
    --predictions data/interim/predictions.parquet \
    --ground-truth data/raw/mibig_ground_truth.tsv \
    --contigs data/interim/analyzed_contigs.txt \
    --output data/processed/benchmark_sharp.json

Reading the output (benchmark.json). docs/METRICS.md is the field-by-field companion — every key in JSON order, written for someone who has not read the code; docs/ARCHITECTURE.md → “Metrics” has the rationale behind the schema:

block what it answers
scope how much of the ground truth was evaluable, and whether scope was explicit or inferred
detection did the tool find the BGC? — min_cluster_frac only
reciprocal the strict symmetric rule, for comparison
nucleotide bp-level agreement; precision says how much extra territory was called
boundary median_prediction_coverage (tightness), split/merge diagnostics

matched_prediction_frac is a lower bound on precision, not precision — the ground truth is incomplete, so an unmatched prediction is unvalidated rather than wrong. There is deliberately no region-level precision or false_positive field.

When a metric is added, renamed, or its meaning changes, update docs/METRICS.md in the same commit — it is the document coworkers read a result with.

Deliberate omissions (do NOT add unless asked)

These were scoped out of the MVP intentionally. Add only when the feature is explicitly needed.

Feature Where it belongs When to add
Evo nucleotide embeddings extract_embeddings.py After ESM-2 baseline is benchmarked
ESM-IF / Foldseek structural embeddings extract_embeddings.py After Evo
GNN embeddings from KG extract_kg_features.py After tabular KG features are validated
Multi-class classification train.py After binary classifier AUROC > 0.85
Ensemble of modality-specific models train.py After multi-class
Asymmetric overlap thresholds metrics.py If team decides one threshold is insufficient
AUROC in benchmark metrics.py Once predict.py scores ALL candidates (not just positives)
Per-class benchmark breakdown evaluate.py When team asks “why is NRPS recall low?”
DeepBGC / antiSMASH per-class breakdown evaluate.py extension When team asks “which BGC class is each tool best at?”
Resumable embedding extraction extract_embeddings.py When datasets exceed ~100k proteins
fp16/bf16 inference model_management.py When running on GPU cluster
BaseStep abstraction new pipeline.py When 3+ steps need to share boilerplate
Logging to file / JSON structured logs config.py When deploying beyond laptop
Docker / Snakemake / Nextflow new When moving to HPC

Key domain facts (don’t get these wrong)

Running the benchmark (current state)

# Smoke test with synthetic data (no real genome needed)
pixi run python scripts/generate_mock_benchmark_data.py \
    --n-clusters 20 --recall-rate 0.7 --n-false-positives 5
pixi run python -m sharp.evaluate \
    --predictions data/mock/predictions.parquet \
    --ground-truth data/mock/ground_truth.tsv \
    --output data/processed/benchmark.json
# → detection recall=0.700 (14/20), matched 14/19 predictions — matches the
#   generator's --recall-rate 0.7 and 5 injected unmatched predictions

# Verify MiBIG 4.0 JSON schema (do once after download)
pixi run python scripts/prepare_mibig_ground_truth.py \
    --inspect data/raw/mibig_json_4.0

# Build real ground truth
pixi run python scripts/prepare_mibig_ground_truth.py \
    --input-dir data/raw/mibig_json_4.0 \
    --output data/raw/mibig_ground_truth.tsv \
    --genus Streptomyces

# BGC Atlas secondary ground truth (noisy — report alongside MiBIG, never alone)
pixi run python scripts/prepare_bgcatlas_ground_truth.py \
    --inspect data/raw/complete-bgcs          # verify schema first
pixi run python scripts/prepare_bgcatlas_ground_truth.py \
    --input-dir data/raw/complete-bgcs \
    --output data/raw/bgcatlas_ground_truth.tsv
#   add --limit N to build against a small subset for dev/tests

# Full competitor comparison (once baseline scripts are written)
# See "Benchmark comparison" section above for full command sequence
This site is open source. Improve this page.