# S(H)ARP S(H)ARP predicts Biosynthetic Gene Clusters (BGCs) in *Streptomyces* and related actinomycetes using **SARP transcription factors as anchors**, combining regulatory context with protein language model embeddings rather than biosynthetic enzyme patterns alone. > GitHub Pages: > GitHub repository: > DAG diagram: Click the links below for the interactive HTML pages: - [DAG diagram (directed acyclic graph)](docs/sharp_dag.html) - Project flowchart, showing inputs-processes-outputs. - [Descriptive page](docs/sharp_pipeline.html) - Practically the same content as the DAG diagram, but presented with a less technical interface. ## Table of Contents - [Quickstart](#quickstart) - [Setting up](#setting-up) - [Pipeline steps](#pipeline-steps) - [Ground truth](#ground-truth) - [Benchmarking](#benchmarking) - [Run once, slice many](#run-once-slice-many--the-benchmarking-approach) - [Benchmarking in a shared clone](#benchmarking-in-a-shared-clone) - [Which genomes, and which slices](#which-genomes-and-which-slices) - [Comparing tools without a ground truth](#comparing-tools-without-a-ground-truth) - [Utilities](#utilities) - [Tests](#tests) - [Directory Structure](#directory-structure) - [Currently Working on](#currently-working-on) ## Quickstart Just cloned the repo? This runs end to end with no data download: ```bash pixi install # set up the environment pixi run pytest # confirm it works # Benchmark against synthetic data — clusters and predictions that overlap # by construction, so the expected number is known in advance. 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 ``` Where to go next: | If you want to... | Read | |---|---| | build the ground truth from MiBiG | [Ground truth](#ground-truth) | | compare against antiSMASH / DeepBGC / GECCO | [Benchmarking](#benchmarking) | | run a S(H)ARP pipeline step | [Pipeline steps](#pipeline-steps) | | move data onto a scratch disk | [Machine-specific settings](#machine-specific-settings-optional-env) | ## Setting up First run: ```bash git clone https://github.com/ecdyzone/sharp cd sharp ``` then install with `pixi` or `conda` ### Using pixi (recommended) Just run: ```bash pixi install ``` ### Using conda Option 1: use `pixi.lock` ```bash conda create --name my-env --file pixi.lock ``` Option 2: use `environment.yml`: ```bash conda env create -f environment.yml ``` ### Machine-specific settings (optional `.env`) Everything works out of the box with no `.env` — data is read from and written to `./data` and the embedding step picks its device automatically. Create one only when this machine differs, e.g. a server that keeps the data on a scratch disk: ```bash cp .env.example .env $EDITOR .env ``` ```bash # .env — gitignored, never committed SHARP_DATA_ROOT=/scratch/$USER/sharp-data # default: /data SHARP_DEVICE=cpu # default: auto (CUDA -> MPS -> CPU) ``` Setting `SHARP_DATA_ROOT` moves `raw/`, `interim/`, `processed/` and `mock/` with it; override one individually with `SHARP_RAW_DIR`, `SHARP_INTERIM_DIR`, `SHARP_PROCESSED_DIR` or `SHARP_MOCK_DIR`. See `.env.example` for all keys. A real environment variable beats the file, so a single run can be redirected without editing anything: ```bash SHARP_DEVICE=cpu pixi run python -m sharp.extract_embeddings ... ``` `.env` holds **machine identity only** — where the data lives, what hardware is present, and credentials once a step needs them. Pipeline parameters (thresholds, model choice, batch sizes) stay on the CLI so a run remains reproducible from the command that produced it. ## Pipeline steps > For reproducing with conda/mamba you have to: > > - `conda activate ` > - run the commands below without `pixi run` ### Extract Embeddings ```bash ## 1. Generate test data pixi run python scripts/generate_mock_data.py --n 100 ## 2. Run the step against it pixi run python -m sharp.extract_embeddings \ --input data/mock/neighborhood_proteins.faa \ --output data/interim/embeddings.parquet ``` ## Ground truth Every benchmark number is measured against one of these two tables. **MiBiG is the primary ground truth** — manually curated, and what any reported result should be based on. **BGC Atlas is secondary and noisy**: its labels are themselves antiSMASH predictions, so agreement with them does not prove correctness. Report it alongside MiBiG, never alone. Build MiBiG first — [Selecting the genome set](#1-selecting-the-genome-set) consumes its output. > **Coverage caveat.** ~53% of MiBiG 4.0 *Streptomyces* entries store no genomic > coordinates, so they cannot be scored and are dropped. The resulting ground truth > is around 400 loci, not 900, and the dropped half skews toward older compound-first > submissions. This affects every tool equally, so it does not bias the comparison — > but absolute recall is "recall over coordinate-resolved MiBiG". See `CLAUDE.md` > for the full analysis. ### Preparing MiBiG Database ```bash # All clusters pixi run python scripts/prepare_mibig_ground_truth.py \ --input-dir data/raw/mibig_json_4.0 \ --output data/raw/mibig_ground_truth.tsv # Or focused on your organism of interest pixi run python scripts/prepare_mibig_ground_truth.py \ --input-dir data/raw/mibig_json_4.0 \ --output data/raw/streptomyces_ground_truth.tsv \ --genus Streptomyces # Bacteria only — drops fungal/plant/animal entries (22% of the # coordinate-resolved set: 327 fungal + 26 plant/animal) pixi run python scripts/prepare_mibig_ground_truth.py \ --input-dir data/raw/mibig_json_4.0 \ --output data/raw/bacterial_ground_truth.tsv \ --exclude-eukaryotes ``` MiBiG's taxonomy block carries only `{name, ncbiTaxId}` — no lineage — so "keep only bacteria" cannot be expressed as a field test. `--exclude-eukaryotes` applies a curated genus deny-list (`EUKARYOTIC_GENERA`, ~100 genera) against the first word of the taxonomy name, and logs which genera it skipped. It is a deny-list rather than an allow-list of bacteria so that an unlisted eukaryote is *kept* (and visible in the log) rather than a novel bacterium vanishing silently. **Not every entry is scoreable.** Of MiBiG 4.0's 3,013 entries, **1,363 (45%) store `location: {from: 0, to: 0}`** — the compound is characterized but the genomic locus is unknown, and a coordinate-based benchmark has nothing to score against. These are dropped on ingest and the count is logged. What remains: | scope | clusters | accessions | |---|---|---| | coordinate-resolved, all genera | 1,634 | 1,420 | | bacteria only (`--exclude-eukaryotes`) | 1,280 | 1,112 | | *Streptomyces* only (`--genus`) | 414 | 352 | Absolute recall figures are therefore "recall over coordinate-resolved MiBiG", not "recall over all known BGCs". The dropped half skews toward older, compound-first submissions — but it is dropped identically for every tool, so the *comparison* stays unbiased. ### Preparing BGC Atlas Database Secondary, noisy ground truth (labels are themselves antiSMASH predictions — report alongside MiBiG, never alone). The dump is 204k antiSMASH `.gbk` files downloaded by `scripts/download_bgc-atlas.sh` (DVC-managed). ```bash # Inspect a few real files first (verify the schema on disk) pixi run python scripts/prepare_bgcatlas_ground_truth.py \ --inspect data/raw/complete-bgcs # Build the TSV (streams over all ~204k files) pixi run python scripts/prepare_bgcatlas_ground_truth.py \ --input-dir data/raw/complete-bgcs \ --output data/raw/bgcatlas_ground_truth.tsv # Develop / test against a small subset without walking 10 GB pixi run python scripts/prepare_bgcatlas_ground_truth.py \ --input-dir data/raw/complete-bgcs \ --output data/interim/bgcatlas_sample.tsv --limit 100 ``` ## Benchmarking Read [How evaluation works](#how-evaluation-works) first — the scope rule it describes governs every command in this section. ### How evaluation works Every tool, including S(H)ARP itself, is scored the same way: its output is converted to a `predictions.parquet` and passed to `sharp.evaluate`, which writes a `benchmark.json`. | block | what it answers | |---|---| | `scope` | how much of the ground truth was evaluable; `explicit` vs `inferred` | | `detection` | *did the tool find the BGC?* (fraction of the cluster covered) | | `reciprocal` | the strict symmetric rule, reported for comparison | | `nucleotide` | bp-level agreement; `precision` = how much extra territory was called | | `boundary` | tightness of the calls, plus split/merge diagnostics | Two things the schema is deliberate about: - **Recall counts only contigs the tool was run on.** Pass `--contigs` (one name per line, or a `.fai`), and give **every tool in a comparison the same file**. Omitted, the scope is inferred from the predictions and a warning is logged — that is optimistic, because a contig analyzed but not called on drops out of the denominator. - **Unmatched is not false.** Ground truth is incomplete by construction, so a prediction with no match is unvalidated rather than wrong. The output reports `matched_prediction_frac` (a lower bound on precision) and `unmatched_prediction_ids` — there is no region-level `precision` field. **[docs/METRICS.md](docs/METRICS.md) explains every field of a `benchmark.json`, in JSON key order** — open a result beside it and read down. Start there if you are interpreting a number rather than producing one. See `docs/ARCHITECTURE.md` → "Metrics — methodological choices" for the full rationale behind the schema. ### Smoke test with mock data Needs no genome and no ground truth — the generator builds clusters and predictions that overlap by construction, so the expected recall is known before the run: ```bash # 1. Generate correlated mock data — clusters and predictions that overlap by construction pixi run python scripts/generate_mock_benchmark_data.py \ --n-clusters 20 --recall-rate 0.7 --n-false-positives 5 # 2. Evaluate pixi run python -m sharp.evaluate \ --predictions data/mock/predictions.parquet \ --ground-truth data/mock/ground_truth.tsv \ --output data/processed/benchmark.json # Output: detection recall=0.700 (14/20), matched 14/19 predictions # — matches the generator's --recall-rate 0.7 and 5 injected extras ``` ### Competitor baselines (antiSMASH / DeepBGC / GECCO) S(H)ARP does **not** run these tools — each has incompatible dependencies and installs into its own isolated pixi env via `scripts/setup_.sh`. You run the tool yourself, then convert its output to `predictions.parquet` and evaluate it exactly like S(H)ARP's own predictions. #### Installing a baseline All three locations the setup scripts write to are set in `.env` (see `.env.example`), so a laptop and a server can differ without editing a tracked file. Defaults apply when `.env` is absent, and each script echoes the paths it resolved: | Key | Default | What | |---|---|---| | `TOOLS_INSTALL_DIR` | `~/.local/src` | one isolated pixi env per tool | | `ANTISMASH_DOWNLOADS_DIR` | `${DATABASES:-~/.local/share}/antismash/databases` | ~10GB reference data | | `DEEPBGC_DOWNLOADS_DIR` | `${DATABASES:-~/.local/share}/deepbgc/data` | ~3GB models + Pfam | The database dirs follow `$DATABASES` when the machine exports it (as the server does) and fall back to `~/.local/share` otherwise, so neither machine normally needs an edit. The `cd ~/.local/src/` commands below assume the default `TOOLS_INSTALL_DIR` — if you override it, use the path the setup script prints when it finishes. #### Running and converting one tool antiSMASH is the worked example; DeepBGC and GECCO differ only where noted below. ```bash # 1. Install a baseline into its own env (~/.local/src//), one-time bash scripts/setup_antismash.sh # 2. Run it yourself, from its own env (or on HPC / in a container) cd ~/.local/src/antismash && pixi run antismash --output-dir # for non-annotated fasta, the code changes a bit: # cd ~/.local/src/antismash && pixi run antismash --output-dir --genefinding-tool prodigal # 3. Convert its output to predictions.parquet (runs in the S(H)ARP env) # (antiSMASH, DeepBGC, and GECCO converters all written) # Inspect first to verify the schema against your actual output: pixi run python scripts/convert_antismash_to_parquet.py --inspect pixi run python scripts/convert_antismash_to_parquet.py \ --input --output data/interim/antismash_predictions.parquet # 4. List the contigs the tool was run on — the denominator of recall. # Build it once from the input genome and reuse it for EVERY tool, or the # recall denominators differ and the numbers stop being comparable. grep '^>' | cut -c2- | cut -d' ' -f1 \ > data/interim/analyzed_contigs.txt # 5. Evaluate against the same ground truth and scope as S(H)ARP pixi run python -m sharp.evaluate \ --predictions data/interim/antismash_predictions.parquet \ --ground-truth data/raw/mibig_ground_truth.tsv \ --contigs data/interim/analyzed_contigs.txt \ --output data/processed/benchmark_antismash.json ``` **DeepBGC** — same shape; S(H)ARP parses `.bgc.tsv` instead: ```bash bash scripts/setup_deepbgc.sh cd ~/.local/src/deepbgc && pixi run deepbgc pipeline --output out pixi run python scripts/convert_deepbgc_to_parquet.py --inspect out pixi run python scripts/convert_deepbgc_to_parquet.py \ --input out --output data/interim/deepbgc_predictions.parquet ``` **GECCO** — same shape, but its `start`/`end` are 1-based inclusive (the one baseline tool needing a coordinate conversion), which the converter applies automatically: ```bash bash scripts/setup_gecco.sh cd ~/.local/src/gecco && pixi run gecco run --genome --output-dir out pixi run python scripts/convert_gecco_to_parquet.py --inspect out pixi run python scripts/convert_gecco_to_parquet.py \ --input out --output data/interim/gecco_predictions.parquet ``` Both then evaluate exactly like step 5 above, with the same `--ground-truth` and the same `--contigs` file, writing `benchmark_deepbgc.json` / `benchmark_gecco.json`. #### Running on Slurm `scripts/run_antismash.sbatch` and `scripts/run_deepbgc.sbatch` wrap the tool invocation with explicit paths and a preflight check. Paths are set at the top of each file — edit them for your machine: ```bash sbatch scripts/run_antismash.sbatch # benchmark genome sbatch scripts/run_antismash.sbatch /path/to/other.fasta sbatch scripts/run_deepbgc.sbatch # benchmark genome sbatch scripts/run_deepbgc.sbatch /path/to/other.fasta squeue -u $USER tail -f deepbgc-scoe-.out ``` Both request CPU only, and no GPU. DeepBGC is Prodigal → `hmmscan` vs Pfam (both CPU-bound, and together the bulk of the runtime) → a small Keras classifier. Only that last stage could use a GPU and it is a rounding error in the total, so a GPU allocation would sit idle — leave the GPU nodes for the ESM-2 embedding step. antiSMASH has no GPU code path at all. DeepBGC's sizing (2 cores / 8G / 2h) comes from `seff` on a real run of this script over `AL645882.2`: **0.86 cores, 1.71 GB, 32m24s**. Despite `hmmscan` dominating the runtime, DeepBGC does not thread it, so the pipeline is serial — an earlier 8-core allocation ran at 10.73% CPU efficiency. antiSMASH's sizing (4 cores / 4G) is measured the same way. It accepts `--cpus` and hands it to its own module scheduler, but that scheduler parallelises very little in practice: `seff` on array job 45315 index 1 — which ran the same `AL645882.2` genome this script defaults to — reported 5.40% CPU efficiency of 16 cores and 1.62 GB peak, i.e. ~1.3 cores of real parallelism. Both `run_antismash.sbatch` and `run_antismash_array.sbatch` are sized from that. The single-genome script keeps a loose 12h walltime because it accepts an arbitrary genome as `$1`, where the array runs a known set. Measure your own runs rather than trusting these numbers on a different genome: ```bash seff sacct -j --format=JobID,State,Elapsed,MaxRSS,ExitCode ``` ### Run once, slice many — the benchmarking approach **This is how every benchmark is produced.** Run each baseline over the broadest genome set you are willing to pay for, *once*, then derive as many benchmark numbers as you like by re-scoping — no recomputation. This works because the expensive step is keyed by genome, not by experiment: ``` ~/projects//out_benchmark// ← one shared pool, keyed by accession │ array tasks skip genomes already done ├─ merge_predictions.py --contigs A ──▶ _predictions_A.parquet ──▶ benchmark_A_.json ├─ merge_predictions.py --contigs B ──▶ _predictions_B.parquet ──▶ benchmark_B_.json └─ merge_predictions.py --contigs C ──▶ ... ``` `--contigs` filters **both** the predictions and the ground-truth denominator (`metrics.py`), so a scope file plus its matching ground truth fully define a benchmark. Producing a new one costs minutes. **Two rules that make this safe:** 1. **Never partition the output pool per experiment.** `OUTROOT` is keyed by accession precisely so two experiments sharing a genome share one directory — that is what makes re-slicing free. Keep the *derived* artifacts separate instead, named per scope: `benchmark_set_/`, `_predictions_.parquet`, `benchmark__.json`. 2. **`--contigs` is mandatory, not optional.** Because the pool holds every genome ever run, `merge_predictions.py` without `--contigs` would sweep in genomes from unrelated experiments. It is what scopes the pool down to one benchmark. A scope is always a **pair** of files from one `select_benchmark_genomes.py --output-dir`: `analyzed_contigs.txt` *and* its `benchmark_ground_truth.tsv`. They must travel together — the ground truth is contig-normalized for that selection (RefSeq/GenBank twins collapsed onto a primary accession), so pairing a scope file with the raw MiBiG ground truth silently drops clusters filed under a twin. ### Which genomes, and which slices The pool we committed to, every scope carved out of it, and the commands that build both, are catalogued in **[docs/BENCHMARK_SCOPES.md](docs/BENCHMARK_SCOPES.md)** — read it before building a scope or submitting an array. In short: the baselines run once over **all 1,280 coordinate-resolved bacterial MiBiG clusters (1,112 accessions)**, BGC-only deposits included. That single run costs ~3 h of antiSMASH and ~1 day of DeepBGC at `%8`, and every scope after it is free: | scope | clusters | what it answers | |---|---|---| | `benchmark_set_strep` | 156 | the headline number | | `benchmark_set_actino` | 575 | where S(H)ARP actually applies — SARPs are an actinomycete family | | `benchmark_set_bact` | ~470-600 | "the entire (bacterial) MiBiG" | | `benchmark_set_pseudomonas` | 61 | negative control — no SARPs, so S(H)ARP should *not* win | | `benchmark_set_smoke` | ~5 | is my setup working? | | per-class slices | NRPS 302, PKS 235, … | which class is each tool best at | The genomes are downloaded once too, into `data/raw/genomes/.fasta`, keyed by accession — so a new scope needs no download of its own either. See also [TODO.md](TODO.md) for the scenarios queued against the shared pool. ### Benchmarking in a shared clone **If someone has already run the baselines, you do not run them again.** On the server the genomes, the output pool and the scope files all live in one checkout, and everyone benchmarks inside it. Producing a number is then a minutes-long merge-and-score, not a multi-day job array. ```bash cd # e.g. ~/projects/igem2026/git-clones/sharp scripts/run_benchmark.sh # no scope: lists the ones that exist scripts/run_benchmark.sh benchmark_set_strep --label $USER ``` That is the whole workflow. `--label` is the only thing you need to remember, and it exists because results are otherwise keyed by scope and tool alone: without it, two people scoring the same scope write the same `benchmark__.json`. With it you get `benchmark__.