Metadata-Version: 2.4
Name: polars-seq-match
Version: 0.1.0
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Developers
Classifier: Programming Language :: Rust
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Programming Language :: Python :: Implementation :: CPython
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Dist: polars>=1.10
License-File: LICENSE
Summary: Fast approximate matching of protein and peptide sequences, as a Polars plugin
Keywords: polars,bioinformatics,protein,peptide,fuzzy-matching
Author-email: Chris Thorpe <drchristhorpe@googlemail.com>
License-Expression: MIT
Requires-Python: >=3.10
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Changelog, https://github.com/drchristhorpe/polars_seq_match/blob/main/CHANGELOG.md
Project-URL: Homepage, https://github.com/drchristhorpe/polars_seq_match
Project-URL: Issues, https://github.com/drchristhorpe/polars_seq_match/issues
Project-URL: Repository, https://github.com/drchristhorpe/polars_seq_match

# polars-seq-match

Fast approximate matching of protein and peptide sequences, as a
[Polars](https://pola.rs) plugin.

Give it two frames — queries and references — and it returns one long frame
telling you which references each query matched and **how close** each match
was: identity, edit distance, and optionally the alignment and a per-residue
list of differences.

The core is domain-neutral. It sees two arrays of strings and returns row
indices plus scores; it knows nothing about alleles, loci or organisms. Domain
rules — which of two equally-good matches wins, which subset of references to
search — stay with the caller, and are expressed as ordinary DataFrame
operations.

```python
import polars as pl
import polars_seq_match as psm

queries = pl.DataFrame({"sample": ["s1", "s2"], "sequence": [seq_a, seq_b]})
references = pl.DataFrame({"allele": ["HLA-A*02:01", ...], "sequence": [...]})

psm.match_frames(queries, references, top_n=3, with_alignment=True)
```

```
┌────────┬──────┬──────────┬──────────┬─────────────┬─────────────┬───────────┐
│ sample ┆ rank ┆ identity ┆ distance ┆ match_type  ┆ mutations   ┆ allele    │
╞════════╪══════╪══════════╪══════════╪═════════════╪═════════════╪═══════════╡
│ s1     ┆ 1    ┆ 1.0      ┆ 0        ┆ exact       ┆             ┆ HLA-A*02… │
│ s2     ┆ 1    ┆ 0.997143 ┆ 2        ┆ approximate ┆ A152E       ┆ HLA-A*02… │
└────────┴──────┴──────────┴──────────┴─────────────┴─────────────┴───────────┘
```

## Install

Requires a Rust toolchain to build from source.

```bash
uv venv --python 3.14
uv pip install maturin polars
maturin develop --release
```

## Two questions, two modes

**`mode="global"`** — how similar are these two whole sequences? This is allele
assignment, protein identification, deduplication.

**`mode="contains"`** — does this short query occur *inside* this reference, and
where? This is epitope mapping and peptide lookup. It reports `match_start` and
`match_end`, and tolerates substitutions.

```python
psm.match_frames(
    pl.DataFrame({"sequence": ["SIINFEKA"]}),
    proteome,
    mode="contains",
    top_n=1,
    with_alignment=True,
)   # -> distance 1, mutations "L280A", match_start/match_end in the protein
```

## Metrics

| `metric` | What it measures | Use it for |
|---|---|---|
| `indel` **(default)** | Edits allowing only insert/delete (a substitution costs 2) | Anything, and specifically anything that must reproduce `Levenshtein.ratio()` |
| `levenshtein` | Classic unit-cost edit distance | When you want "one substitution = distance 1" |
| `hamming` | Positional mismatches, equal lengths only | Fixed-length peptide sets (9-mers) |
| `jaro_winkler` | Prefix-weighted similarity | Name-like strings, not really sequences |

`indel` identity is *numerically identical* to Python `Levenshtein.ratio()` and
`rapidfuzz.fuzz.ratio() / 100`. That is tested against both packages, so
migrating an existing pipeline does not silently move your numbers.

## Reading the output

Always present:

| Column | Meaning |
|---|---|
| `identity` | Normalised similarity in `[0, 1]` |
| `distance` | Edit operations (null for `jaro_winkler`) |
| `match_type` | `exact` or `approximate` |
| `rank` | 1-based, best first |
| `group_size` | How many reference rows share this exact sequence |
| `exhaustive` | False if a truncating prefilter may have missed something |

With `with_alignment=True` you also get `query_aligned`, `reference_aligned`,
`n_substitutions` / `n_insertions` / `n_deletions`, and `mutations` — a compact
`A152E,D260del,97insK` summary in 1-based reference coordinates.
`psm.parse_mutations(result)` explodes that into one row per difference.

Queries that match nothing still appear, once, with null match columns. The
result reads as a left join.

## Identical reference sequences

Reference rows sharing a sequence are grouped. `top_n` counts **distinct
sequences**, and by default each contributes one row — the earliest in
reference order — with `group_size` telling you how many names it stands for.
That default exists because one MHC protein can carry hundreds of allele names.

```python
result["group_size"][0]            # 40 names share this protein
index.group_frame(reference_index) # all 40 rows, on demand
psm.match_frames(..., expand_groups=True)   # or fan them out up front
```

### Ties resolve by reference row order

The library has no opinion about which of two equally-good matches wins — it
takes the one appearing **earlier in your reference frame**. Sort the reference
frame by your own priority first and the matcher honours it:

```python
references = references.sort(by=my_priority)   # e.g. allele number ascending
```

That single rule is how `histo_match`'s "lowest allele number wins" is
reproduced without the library ever learning what an allele number is.

## Matching many batches

Building the index — deduplicating, ordering by length, building trimer
postings — is the expensive part. Do it once:

```python
with psm.SequenceIndex(references, kmer_size=3) as index:
    a = index.match(batch_one, top_n=5)
    b = index.match(batch_two, top_n=5)
```

Metric and mode are query-time options, so one index serves all of them.

## Speed, and what to reach for

Measured on 14,200 references × 1,000 queries of ~350 residues, `top_n=5`
(`python benchmarks/bench.py`):

| | Time | Exact? |
|---|---|---|
| Python loop over `Levenshtein.ratio` | 25.7 s | yes |
| `rapidfuzz.process.cdist` (C++, multithreaded) | 2.9 s | yes |
| this plugin, default | 2.8 s | yes |
| this plugin, trimer prefilter + `max_distance=20` | 0.9 s | **yes, provably** |
| this plugin, trimer prefilter, top-512 seeds | 0.7 s | no (flagged) |

The default path is **exhaustive** — it returns the true top-N, and that is
checked against brute force in the test suite. It gets its speed from an exact
length bound, which prunes hard when references vary in length (the same run
drops to 0.66 s versus `cdist`'s 1.56 s) and not at all when they are all the
same length — the realistic case for full-length MHC alleles.

For that flat-length case, reach for the trimer prefilter *with a distance
bound*, which stays provably exact:

```python
index = psm.SequenceIndex(references, build_kmer_index=True)
index.match(queries, prefilter="kmer", max_distance=20,
            prefilter_candidates=len(references))
```

Without a `max_distance` the prefilter is a heuristic and can miss a better
match. It says so: the `exhaustive` column comes back false. Nothing in this
library approximates silently.

Synthetic data is easy to make flattering, so the same run against a real MHC
class I reference set is worth reporting — 25,825 sequences over 13,748 distinct
proteins of 248-275 residues, 1,000 queries each carrying three substitutions.
(`benchmarks/bench.py --real` reproduces this against a local copy of that data;
it is third-party and not distributed here, and the flag exits cleanly without
it.)

| | Time | Queries/s | Exact? |
|---|---|---|---|
| exhaustive | 1.75 s | 572 | yes |
| exhaustive + `max_distance=30` | 1.67 s | 599 | yes |
| trimer prefilter + `max_distance=30` | 1.05 s | 950 | yes |

Every row there is exhaustive: the distance bound keeps the trimer prefilter
provably lossless, so the 1.7x is free. Index build is a further 0.12 s, paid
once.

[EXPLANATION.md](https://github.com/drchristhorpe/polars_seq_match/blob/main/EXPLANATION.md) covers why all of this works.

## Development

```bash
uv sync                                  # builds the extension and installs dev deps
uv run pytest                            # 62 tests
uv run python benchmarks/bench.py

cargo fmt --check
cargo clippy --all-targets -- -D warnings

# Rust unit tests link against libpython, so it needs to be on the loader path:
LD_LIBRARY_PATH=$(python -c 'import sysconfig;print(sysconfig.get_config_var("LIBDIR"))') \
  PYO3_PYTHON=$(which python) cargo test --release
```

The tests that matter most are in `tests/test_search_correctness.py`: they check
that every optimisation returns exactly what a brute-force scan would.

CI runs the Python suite on 3.10 and 3.14, the Rust unit tests, rustfmt and
clippy, and installs the built wheel on the abi3 floor to check the packaged
artefact rather than only the source tree. Releases are tag-driven — see
[RELEASING.md](https://github.com/drchristhorpe/polars_seq_match/blob/main/RELEASING.md).

## Limitations

- No substitution matrices (BLOSUM/PAM) or affine gaps. For divergent
  homologues, edit distance is the wrong model — use a real aligner.
- `jaro_winkler` has no pruning bound and always scans everything.
- Sequences must be ASCII. Non-ASCII input is rejected rather than compared
  byte-wise and quietly wrong.
- `kmer_size` is 2–4 (the posting index is dense).

## See also

- [PLAN.md](https://github.com/drchristhorpe/polars_seq_match/blob/main/PLAN.md) — design rationale and what changed while building it
- [EXPLANATION.md](https://github.com/drchristhorpe/polars_seq_match/blob/main/EXPLANATION.md) — how sequence comparison works, and how it
  maps onto Polars
- [CHANGELOG.md](https://github.com/drchristhorpe/polars_seq_match/blob/main/CHANGELOG.md)

