Metadata-Version: 2.5
Name: stdb-mlst
Version: 0.1.2
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` — mmseqs2
- `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
```

- 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), `--blast-bin` (defaults to the
  `STDB_BLAST_BIN` env var or `/home/shichang/.local/micromamba/envs/etoki/bin`).
- 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.

## 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
```
