Metadata-Version: 2.4
Name: prismalign
Version: 0.2.11
Summary: N-color (2-color / 3-color / 3-nt) nucleotide-conversion alignment engine with pluggable backends
Author-email: Chang Ye <yech1990@gmail.com>
License-Expression: GPL-3.0-only
Keywords: bioinformatics,nucleotide conversion,bisulfite,SLAM-seq,m6A,alignment,epigenetics
Classifier: Programming Language :: Python :: 3
Classifier: Operating System :: POSIX :: Linux
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: bwamem>=0.0.56
Requires-Dist: pysam>=0.22
Provides-Extra: mappy
Requires-Dist: mappy>=2.24; extra == "mappy"
Provides-Extra: minibwa
Requires-Dist: minibwa>=0.1.7; extra == "minibwa"
Provides-Extra: bwamem2
Requires-Dist: bwamem2>=0.1.0; extra == "bwamem2"
Dynamic: license-file

# Prismalign

**N-color nucleotide-conversion alignment engine** with pluggable backends.

Prismalign maps sequencing reads from any nucleotide-conversion chemistry
(bisulfite-seq `C→T`, SLAM-seq `T→C`, m6A / A-to-I `A→G`, MK/KM dual-base,
or a custom 3rd channel) using a **HISAT-3N-style** strategy:

1. build a *converted* reference index (`scheme.ref_from → ref_to`)
2. transform each read per color channel and align it to the converted index
   via a pluggable backend (**bwamem** by default; **WFA2-lib**;
   **minimap2/mappy**; **minibwa**, all optional)
3. **re-score every hit against the *original* reference** so that real
   conversions are rewarded (not counted as mismatches), emitting a
   color-correct `MD` plus per-channel `Y`/`Z` counts in BAM tags.

All per-read heavy kernels are native C (BWA-MEM / WFA2-lib / minimap2 /
minibwa); the Python layer is a thin, friendly wrapper.

## Install

```bash
pip install -e .              # bwamem + built-in WFA2 C backends
pip install -e "./[mappy]"    # + minimap2 backend
```

## Usage — Python (clean wrapper)

```python
import prismalign as ps

# one-shot mapping -> BAM (builds indexes, maps, cleans up)
ps.map_reads("reads.fq", "ref.fa", "out.bam",
             scheme="MK", backend="bwamem", threads=4)

# object API / reuse
with ps.NColorMapper(scheme=ps.BS, backend="bwamem") as mapper:
    mapper.map_file("reads.fq", ref_files=["ref.fa"], output_files=["bs.bam"])
```

## Usage — CLI

```bash
# classic two-color (MK: A->G + C->T); --backend auto picks the fastest
# importable backend (minibwa > mappy > bwamem)
prismalign map -s MK -r ref.fa -o out.bam reads.fq

# force the fastest native backend (bwa-mem2 speed tier)
prismalign map -s MK --backend minibwa -r ref.fa -o out.bam reads.fq

# bisulfite-seq (3-nt single channel C->T)
prismalign map -s BS -r genome.fa -o bs.bam --index-dir idx reads.fq

# parallel (2 copies of the reads, byte-identical output to -t 1)
prismalign map -s MK -r ref.fa -o out.bam -t 4 reads.fq

# drive an external CLI tool via the fast batched path (--adapter, never --backend)
prismalign map -s MK --adapter bwa-mem2 -r ref.fa -o out.bam reads.fq

# list built-in schemes
prismalign schemes
```

## Schemes

| name  | reference index | channels | use case |
|-------|-----------------|----------|----------|
| `MK`  | `AC→GT`         | 2        | dual-base conversion A→G + C→T (classic two-color) |
| `KM`  | `GT→AC`         | 2        | reverse of MK |
| `BS`  | `C→T`           | 1        | bisulfite-seq (3-nt) |
| `SLAM`| `T→C`           | 1        | SLAM-seq |
| `A2G` | `A→G`           | 1        | m6A / A-to-I editing |
| `THREE`| `AC→GT`        | 3        | three-color demo (add your 3rd base pair in `schemes.py`) |

## Python API

```python
from prismalign import NColorMapper, BS

mapper = NColorMapper(scheme=BS, backend="bwamem", index_dir="idx")
mapper.map_file(r1_file="reads.fq", ref_files=["genome.fa"],
                output_files=["out.bam"])
```

## Backends vs Adapters — two integration layers

Prismalign plugs in aligners at **two layers**:

1. **Backends** — **in-process** (compiled C / a Python binding). The engine's per-read
   `map_one` calls `backend.align(seq)` directly. Selectable via `--backend`.
2. **Adapters** — **subprocess** wrappers for external *command-line* aligners. They
   map reads at three granularities (`map_read`/`map_batch`/`map_file`), unified on the
   `CliAdapter` base; *not* `--backend`-selectable.

### Backends (in-process, `--backend`)
| backend | engine | notes |
|---------|--------|-------|
| `bwamem` | BWA-MEM via the `bwamem` package | default, fast C backend (SE + PE) |
| `minibwa` | **lh3/minibwa** (bwa-mem successor) via **PyO3 pip binding `minibwa`** (fg-labs) | ~2-3x faster than bwa-mem; `pip install minibwa` (SE + PE) |
| `mappy` | minimap2 via `mappy` | official minimap2 Python binding (SE + PE) |
| `bwamem2` | **bwa-mem2** native in-process via the `bwamem2` Cython binding | correct, but **per-read**; `pip install 'prismalign[bwamem2]'` (standard, non-free-threaded CPython) |
| `wfa2` | **WFA2-lib** (vendored v2.3.6, MIT) compiled in-process | exact gapped (indel-aware) wavefront alignment; SE only |

### Adapters (subprocess, `prismalign.adapters`)
All built on the shared `CliAdapter` base, which maps reads at three granularities:

| method | granularity | notes |
|---|---|---|
| `map_read(seq)` | 1 read | compat / debug (slow: one subprocess per read) |
| `map_batch(reads)` | 1 batch (`(name,seq,qual)` list) | one tool invocation |
| `map_file(fastq, batch_size, threads, workers)` | whole FASTQ (plain/.gz) | **path prismalign reads + batches itself** (default `batch_size=8192` = 8×bwa-mem2's internal 512); `workers>1` runs batches **concurrently** (each an independent tool subprocess — the real throughput lever, since a tool's own threads saturate on a shared index) |

| adapter | tool | index |
|---|---|---|
| `BwaMemAdapter` | `bwa mem` | `bwa index` |
| `BwaMem2Adapter` | `bwa-mem2 mem` | `bwa-mem2 index` |
| `Bowtie2Adapter` | `bowtie2` | `bowtie2-build` |
| `Minimap2Adapter` | `minimap2 -a` | none (reads FASTA directly) |
| `Hisat2Adapter` | `hisat2` | `hisat2-build` |
| `StrobealignAdapter` | `strobealign` | `.sti` |
| `SamAdapter` | any SAM mapper | (generic; you supply the command template) |

Each adapter needs its binary on PATH (or an env var: `BWA_BIN`, `BWA_MEM2_BIN`,
`BOWTIE2_BIN`, `MINIMAP2_BIN`, `HISAT2_BIN`, `STROBEALIGN_BIN`).

**Throughput tip:** for a large job use ``map_file(path, batch_size=8192,
workers=cpu/2..cpu, threads=1)`` — batching at 8192 (a multiple of the tool's
internal 512) amortizes the subprocess startup, and ``workers`` parallelizes
batches across **separate tool processes** (each with its own memory bandwidth;
the tool's own ``-t`` saturates, so parallelize via processes). `.gz` input is
read with an internal fast library (`isal`/`xopen`/`rapidgzip`) when available.

### bwa-mem2 appears at BOTH layers (intentional)
- **`BwaMem2Backend`** (backend) = bwa-mem2 **in-process**, per-read — for embedding / `--backend bwamem2`.
- **`BwaMem2Adapter`** (adapter) = bwa-mem2 **batched** CLI (`map_file`) — for throughput.

For large references, bwa-mem2's ~2-3x (from SIMD FM-index search) shows up; build it
**cleanly for AVX2/AVX-512** (stale object files cause SIGILL). The in-process backend
is per-read (batch the calls, or use the adapter, for throughput).

Full inventory — including where each wrapper lives — in [`docs/backends.md`](docs/backends.md).

## Speed & IO

* **Auto backend**: `--backend auto` (explicit) picks the fastest *importable*
  native backend — `minibwa` (~2-3x BWA-MEM, the bwa-mem2 speed tier), else
  `mappy` (in-process), else `bwamem`. (Default is `bwamem`.) The chosen
  backend is printed at startup (`[prismalign] backend auto=minibwa`).
* **Process parallelism is the lever**: `-t/--threads N` maps reads in an
  ordered fork+COW process pool (any backend); batches are drained in read
  order so the BAM is **byte-identical** to `threads=1`. Processes (each with
  its own index) scale; bwa-mem2's internal `-t` threads do *not* (they share
  one index and are memory-bound). So scale **processes/workers**, not
  bwa-mem2 threads.
* **Reduced repeated IO**: references are copy+converted **once** even when
  reused across layers (cache keyed by path+scheme); per-hit reference fetch
  is cached in memory for small contigs (RNA/transcript references), so only
  one indexed read per contig.

> **Throughput ceiling.** The heavy per-read work (BWA-MEM / minibwa / minimap2
> / WFA2 kernels) is native C. With BWA-MEM the runtime is essentially that
> kernel's throughput — prismalign's glue (conversion, re-scoring, tag
> emission) adds only a small fraction. For the largest runs, use `minibwa`
> (runs on standard CPython; no free-threaded-3.14 wheel) or `mappy`; both are
> far faster than BWA-MEM.

## Limitations (v0.2.x)

* paired-end is supported natively by the `bwamem`, `mappy` and `minibwa`
  backends, and by the `bwa-mem2` CLI adapter (batched `map_batch` PE);
  `wfa2` and the remaining subprocess adapters (`sam`/`strobealign`/
  `bowtie2`/`minimap2`/`bwa`) are single-end.
* hierarchical (layered) mapping uses the `PLAIN` identity scheme for
  non-converted short-RNA references. `minibwa` is the fastest native backend
  but needs standard (non-free-threaded) CPython ≤ 3.13, so on 3.14t `mappy`
  is the fastest available.
