Metadata-Version: 2.5
Name: stdb-mlst
Version: 0.1.8
Summary: Local MLST/cgMLST ST database manager: SQLite-backed CLI for importing allele sequences and ST profiles
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.11
Provides-Extra: dev
Requires-Dist: pytest>=8.0.0; extra == 'dev'
Requires-Dist: ruff>=0.5.0; extra == 'dev'
Provides-Extra: gp
Requires-Dist: numpy; extra == 'gp'
Requires-Dist: scikit-learn; extra == 'gp'
Description-Content-Type: text/markdown

# stdb

A local MLST/cgMLST ST database manager: pure CLI built on the Python standard
library `sqlite3` — no network dependency, no server process.

Each species/version's ST library is **one folder + one SQLite file**:

```
data/
└── vparahaemolyticus_3/
    └── library.db          # all library data (alleles + profiles + meta + import log)
```

## Installation

### Conda environment (recommended)

[`environment.yaml`](environment.yaml) provisions a `stdb` environment with
the package and all third-party tools it uses — blast
(blastn/makeblastdb/blastp for `type-genome`, including the protein recovery
layer), mmseqs2 (`init --cluster` / `cluster-refs`) and numpy/scikit-learn
(the exact EToKi GP in `select-loci`):

```bash
micromamba env create -f environment.yaml   # or: conda env create -f environment.yaml
micromamba activate stdb
```

The third-party binaries are installed directly on `PATH`, which is the final
step of the `--blast-bin` / `--mmseqs-bin` resolution chains — no extra
configuration needed.

### PyPI

The core package is pure stdlib with zero runtime dependencies:

```bash
pip install stdb-mlst
```

Tool requirements per feature (binaries via `PATH`, or explicit
`--blast-bin` / `--mmseqs-bin`):

- `type-genome` — blastn/makeblastdb/blastp
- `init --cluster` / `cluster-refs` / `homologs` / `build` / `seed-loci`
  (multi-reference mode) — mmseqs2
- `pangenome` (experimental) — PEPPAN (`pip install PEPPAN`)
- `select-loci` — scikit-learn optional (`pip install numpy scikit-learn` or
  `pip install 'stdb-mlst[gp]'`); without it, variability control falls back
  to a linear-regression proxy

### Development

```bash
pip install -e '.[dev]'     # package + pytest + ruff
```

## Data directory

The data root is resolved in this order:

1. CLI flag `--data-dir DIR`
2. Environment variable `STDB_DATA_DIR`
3. Default `./data`

## Core features

### 1. Initial import (allele sequences + allele profiles)

```bash
stdb init vparahaemolyticus_3 \
  --alleles test/allele.fna \
  --profile test/alleles.st.profile.gz \
  --type cgmlst --species "Vibrio parahaemolyticus" --data-dir data
```

- Idempotent: re-running never duplicates data (exact sequence dedup).
- Allele numbers keep the official numbers (digits in the FASTA headers, e.g.
  `>VP0001_7` → allele_id=7).
- Profiles support gzip; both `#Genome`/`#ST` (cgMLST) and `ST` (pubMLST
  7-gene) headers are auto-detected.
- Missing-allele markers (`-`, `~`, empty, ...) are normalized to 0.
- Same number + different sequence conflicts are counted and skipped without
  blocking the import.

### 2. New allele sequence import (dedup)

```bash
stdb add-alleles vparahaemolyticus_3 new_alleles.fna --data-dir data
```

- Exact dedup by `(locus, sequence)`: existing sequences are skipped.
- New sequences get the locus's current max + 1 in file order (header numbers
  are not trusted and ignored).

### 3. New allele profile import (new ST assignment)

```bash
stdb add-profiles vparahaemolyticus_3 novel_profiles.tsv --data-dir data
```

- Rows with an official ST number (e.g., newly released pubMLST STs) are stored
  under that number; same number + different profile is an error with full
  rollback.
- The same profile under two different official ST numbers (pubMLST historical
  legacy data) → the first-stored ST wins, the later one is recorded as
  `aliased`, skipped, and reported — the import is not blocked.
- Rows without an ST number → exact match against existing STs; on a hit the
  ST is returned, otherwise `max(ST)+1` is assigned (recorded in `new_sts`).

### 4. Genome import (returns allele profile and ST number)

```bash
stdb add-genome vparahaemolyticus_3 genome_alleles.fna --data-dir data
```

- Input: one genome's allele sequences (one sequence per locus; omit missing
  loci).
- Sequences are matched against existing alleles by `(locus, sequence)`, new
  sequences get auto-assigned numbers (header numbers not trusted).
- Returns the genome's allele profile (missing loci = 0) and ST number:
  - profile hits an existing ST → the original number (`st_new: false`);
  - no hit → `max(ST)+1` stored (`st_new: true`).
- Sequence and profile changes land in a single transaction.

### 5. Genome typing (reference-allele based, blastn; report-only by default)

```bash
# report only (no DB writes)
stdb type-genome vparahaemolyticus_3 assembly.fna --data-dir data
# commit after review (new alleles + profile + ST stored)
stdb type-genome vparahaemolyticus_3 assembly.fna --commit --data-dir data
# several genomes at once, typed concurrently, one ST per genome
stdb type-genome vparahaemolyticus_3 g1.fna g2.fna g3.fna --jobs 8 --data-dir data
```

Multi-genome runs execute the per-genome pipelines concurrently (`--jobs`,
default 4; total cores ≈ jobs × threads; concurrent `--commit` writes are
serialized by the library write lock) and report one compact record per
genome in input order — `{genome, st, st_new, calls}`. Single-genome runs
keep the original detailed per-locus output.

- Input: genome assembly FASTA (gzip supported). Library reference alleles are
  the query, the genome is the blastn database (EToKi MLSType approach and
  tuned parameters).
- Per locus, the best hit is taken and the genomic interval extracted:
  - sequence exactly matches an existing allele → `matched`;
  - hit with identity ≥ `--min-identity` (default 0.9) → `novel`
    (committable new allele);
  - hit below threshold → `low_identity` (report-only, never committed);
  - no hit → `missing` (0 in the profile).
- Returns the allele profile and ST (predicted number when no ST matches).
- Only `--commit` writes: new alleles numbered max+1, new profile/ST stored in
  one transaction (write lock + conflict retry); `low_identity`/`missing` are
  never committed.
- Protein recovery layer (on by default, `--no-protein` disables; the blastp
  equivalent of EToKi's dualBlast protein channel): loci blastn left as
  `missing`/`low_identity` are recovered via blastp protein search (reference
  alleles and genome ORFs are each translated in 6 frames; frame/strand-aware
  nucleotide coordinate mapping + EToKi's lookForORF start/stop codon gate).
  Recovered classification:
  - extracted interval (after ORF trimming) md5-exact → `matched`;
  - protein identity ≥ `--min-identity` → `novel` (committable);
  - incomplete ORF (frameshift/truncation, no start or stop codon) → `pseudo`
    (report-only, never committed — mirrors EToKi's refusal of low-quality
    alleles);
  - otherwise → `low_identity`.
  The protein layer builds reference proteins only for blastn-unresolved loci
  (O(unresolved loci) memory).
- Other options: `--min-coverage` (default 0.8), `--loci VP0001,VP0003`
  (subset typing), `--threads` (default 8), `--jobs` (concurrent genomes,
  default 4), `--blast-bin` (defaults to the `STDB_BLAST_BIN` env var,
  `~/.local/micromamba/envs/etoki/bin` when present, or `PATH`).
- Reference sequences default to `<library>/type_refs.fna` (mmseqs cluster
  representatives, below); falls back to a full alleles export when the file
  is absent; `--full-refs` forces the full export.

### 6. Typing reference clustering (mmseqs; optional, speeds up typing)

```bash
# cluster during init (--cluster takes an optional identity, default 0.96)
stdb init vparahaemolyticus_3 --alleles alleles.fna --profile alleles.st.profile.gz --cluster
# rebuild the reference separately (e.g., after adding many alleles)
stdb cluster-refs vparahaemolyticus_3 alleles.fna [--identity 0.96]
```

- One global mmseqs easy-cluster run over the input alleles.fna (linclust
  handles millions of sequences; per-locus invocations are infeasible on
  loaded machines, and cross-locus misclusters do not affect typing
  correctness — hits are still resolved by md5 exact matching).
- Each cluster's representative is mapped back to a real DB allele_id via md5
  point queries (representatives not in the library are dropped), written
  atomically to `data/<lib>/type_refs.fna` (O(1) memory).
- Typing's blastn query count drops from all alleles to the representative
  count, significantly speeding up `type-genome`; a lower `--identity` means
  fewer representatives and faster typing (default 0.96 preserves
  discriminative power).
- mmseqs lookup: `--mmseqs-bin` > `STDB_MMSEQS_BIN` > PATH > etoki env.

### 7. cgMLST allele selection (select-loci; mirrors EToKi cgMLST.py)

```bash
# library mode: select against the library's profiles/alleles, write scheme files
stdb select-loci vparahaemolyticus_3 [--output prefix] \
    [--genepresence 0.95] [--intactcds 0.94] [--genelength 0] [--oddratio 3.0]
# file mode: select directly on input files (no library needed)
stdb select-loci --alleles alleles.fna --profile profiles.tsv [--output prefix]
# select first, then import — one step
stdb init vp_selected --alleles alleles.fna --profile profiles.tsv --select
```

- Per-allele quality status (EToKi semantics): 2=pseudo (length not a
  multiple of 3), 3=internal stop codon, 4=no start codon, 5=intact CDS,
  6=with stop codon.
- Three filter rounds (same structure and defaults as EToKi): loose round
  (presence ≥0.5 / intact ≥0.5, drop STs with genome presence <0.4) →
  medium round (0.8/0.8, ST threshold 0.6) → strict round (default presence
  ≥0.95 / intact ≥0.94 / average length / variability control).
- Variability control: EToKi uses an sklearn GP regression (RBF+WhiteKernel
  kernel, 9 restarts) to drop loci whose variation is significantly above the
  length expectation; stdb uses the **same GP** when sklearn is available
  (`pip install scikit-learn` into the venv; the pExp/#Sigma/CI columns match
  EToKi), and degrades to a linear-regression proxy with a stderr warning
  otherwise. GP fitting is slow (several minutes on the real library; EToKi's
  own comments note "can take a long time").
- Outputs (`<prefix>.alleles.fna` / `<prefix>.profiles.tsv` /
  `<prefix>.cgMLST`): alleles.fna holds all sequences of the kept loci;
  profiles.tsv holds the surviving STs restricted to the kept loci (dedup:
  identical profiles keep the smallest ST number) — both files can be fed
  directly into a new `init`. `--select` runs "select first, then import"
  inside init.
- Library mode is read-only; status precomputation over the 1.38M alleles
  takes a few dozen seconds.

### 8. Pan-genome loci discovery (experimental; wraps PEPPAN)

```bash
# homolog analysis over annotated genomes -> gene families become loci
stdb pangenome --output pan genome1.gff.gz genome2.gff.gz ... \
    [--orthology nj|ml|sbh] [--match-identity 0.65] [--min-cds 120] \
    [--clust-identity 0.9] [--clust-match-prop 0.8] [--threads 8]
# turn the discovered loci into a library
stdb init species_pan --alleles pan.allele.fna --profile pan.profiles.tsv
```

- Homologous-gene analysis is delegated to the external
  [PEPPAN](https://github.com/zheminzhou/PEPPAN) tool (`pip install PEPPAN`;
  needs mmseqs2/diamond/blastn at runtime, bundled with its pip package). stdb
  only orchestrates and converts: each orthologous gene family becomes one
  **locus** of a new scheme.
- PEPPAN resolution chain: `--peppan-bin` > `STDB_PEPPAN_BIN` > PATH.
- Outputs at `<prefix>.*`:
  - `.allele.fna` — PEPPAN's pan-genome allelic sequences; headers are
    `family_allele`, already stdb's `locus_allele` convention;
  - `.profiles.tsv` — binary presence/absence profile (one ST row per genome,
    1/0 per locus);
  - `.genomes.tsv` — ST ↔ genome-name mapping sidecar;
  - `.PEPPAN.gff` — PEPPAN's full re-annotation, kept verbatim.
- Semantics: pseudogene copies still mark presence; merged predictions
  (`inference=ortholog_group:a,b`) mark every listed family; removed/untrusted
  predictions (misc_feature lines) are ignored.
- PEPPAN runs inside a temporary work directory; only the four artifacts above
  land next to the output prefix.

### 9. Building a cgMLST scheme from scratch (experimental; seeded workflow)

```bash
# 1. seed loci from annotated reference genome(s) — several references pool
#    their genes into one union locus set (shared genes merge across
#    references, their sequence variants become alleles; multicopy genes and
#    genes < 120 bp are pre-filtered away)
stdb seed-loci ref1.gff.gz ref2.gff.gz --output seeds   # single GFF: + --fasta genome.fna
# 2. create the library with the seeds (no profile yet)
stdb init species_build --alleles seeds.fna --data-dir data
# 3. extract each genome's alleles (blastn + protein recovery layer)
for g in genomes/*.fna; do stdb type-genome species_build "$g" --commit --data-dir data; done
# 4. flag cross-locus homologous (paralog-suspect) loci
stdb homologs species_build --data-dir data         # -> species_build.homologs.tsv
# 5. statistical optimization of the scheme
stdb select-loci species_build --output selected --data-dir data
# 6. rebuild a clean library from the optimized scheme
stdb init species_cgmlst --alleles selected.alleles.fna --profile selected.profiles.tsv --data-dir data
```

- `seed-loci` takes every CDS of the annotated reference genome(s) as
  candidate loci and writes `<prefix>.fna`. Multicopy seeds (identical
  sequence occurring twice or more within one genome) are dropped entirely —
  typing them would hit several genomic sites and flip between copies
  arbitrarily; short seeds (< `--min-len`, default 120 bp) are dropped.
  Sequences come from the GFF's embedded `##FASTA` section, `<gff>,<fasta>`
  pairs, or (single GFF) `--fasta`.
- Multi-reference mode clusters the pooled candidates (`--identity`, default
  0.95) so the same gene from different genomes collapses into one locus whose
  distinct copies become alleles 1..n; locus names follow the first reference
  carrying the gene. Needs mmseqs2; single-reference mode is pure stdlib.
- `--locus-prefix GAS_` renames loci to systematic ids (`GAS_00001..`) in
  deterministic input order instead of keeping the annotation names; the
  old->new mapping is written to `<prefix>.names.tsv`.
- `homologs` clusters every library allele globally with mmseqs (default
  identity 0.9 / coverage 0.8 on both sequences) and reports clusters whose
  members span more than one locus to `<lib>.homologs.tsv` — paralog suspects
  whose copies cross-hit each other during typing. Review/drop those loci
  before `select-loci`.
- Step order matters: homolog filtering before `select-loci` keeps
  statistically attractive but mutually interfering loci out of the scheme;
  step 3 is embarrassingly parallel across genomes (writes are lock-serialized).

### 10. One-command database construction (experimental; `build`)

```bash
stdb build species_cgmlst \
    --refs ref1.gff.gz ref2.gff.gz \
    --genomes genomes/*.fna \
    [--locus-prefix GAS_] [--seed-identity 0.95] [--homolog-identity 0.9] \
    [--skip-homologs] [--skip-select] [--cluster] \
    [--genepresence 0.95] [--intactcds 0.94] [--oddratio 3.0] \
    [--threads 8]
```

Runs the whole seeded workflow of section 9 in one step:

1. multi-reference seeding over `--refs` (union locus set, sequence variants
   become alleles);
2. a working library `<name>.work` is initialized and every `--genomes` FASTA
   is typed into it (`--commit`);
3. `homologs` flags cross-locus paralog-suspect families and they are excluded
   (`--skip-homologs` disables);
4. `select-loci` optimizes the scheme statistically (`--skip-select` disables);
5. the final library `<name>` is rebuilt from the optimized scheme and every
   `--genomes` FASTA is re-typed against it under its own conditions — the
   stored STs are reproducible by construction when genomes are re-typed
   later; `--cluster` pre-builds its typing refs like `init --cluster`.

Artifacts: `<name>` is the deliverable; `<name>.work` keeps the full
unfiltered data together with `homologs.*` (what was excluded) and
`selected.*` (the intermediate optimized scheme files). Both target names must
not exist yet.

## Read-only queries

```bash
stdb list                              # list all libraries
stdb info <lib>                       # library info (counts, scheme, recent imports)
stdb stats <lib> [--top 10]           # per-locus stats (allele counts, etc.)
stdb profile <lib> <ST>             # profile by ST number
stdb find <lib> --locus VP0001=7 --locus VP0003=30   # exact ST lookup by allele list
stdb seq <lib> --locus VP0001 --allele 7             # sequence by locus+number
```

All commands print JSON to stdout; errors go to stderr.

## Numbering rules

| Scenario | Rule |
| --- | --- |
| Initial allele import | allele_id = FASTA header number (official) |
| New allele import | allele_id = locus current max + 1 |
| Initial profile import | official ST numbers from the profile are kept |
| New profile (no ST number) | reuse the matching ST, else ST = max + 1 |

The upper bound of official ST numbers is recorded in meta as
`max_official_st` after init and subsequent official-ST imports (visible via
`stdb info`); locally added ST numbers are always above it, and
`local_st_count` is the number of local novel STs.

Each locus's max allele number is maintained in the `allele_max_ids` table
(auto-updated in the same transaction as imports, auto-calibrated against the
real data before imports, directly queryable).

## Concurrency safety

Parallel `stdb` invocations (multiple processes/threads writing the same
library) are safe:

- **WAL mode**: readers never block writers.
- **Write-lock serialization**: writes to the same library (`init` /
  `add-alleles` / `add-profiles`) queue on the file lock
  `<lib-dir>/write.lock`; later calls wait automatically — no numbering race
  (two calls computing the same `max+1`) and no `database is locked` failures.
- **Conflict fallback**: even if the file lock is bypassed (e.g., another tool
  writing the DB directly), transaction conflicts roll back and retry
  (max 3, re-reading DB state each attempt), converging on
  duplicate/matched for concurrently written data.
- Example: parallel batches of new sequences/profiles execute one after
  another, producing exactly the serial result.

## Tests

```bash
python -m pytest tests/ -v          # synthetic-data unit tests
python -m ruff check src tests      # lint
```

Optional real-data smoke test (test/ holds 1.6GB of Vibrio parahaemolyticus
cgMLST data):

```bash
stdb init vparahaemolyticus_3 --alleles test/allele.fna --profile test/alleles.st.profile.gz --data-dir data
```

`tests/test_real_genomes.py` additionally runs the seeded-pipeline smoke test
(seed-loci -> init -> type-genome -> homologs) on real GAS genomes when the
directory pointed to by `STDB_REAL_GENOMES` exists.
