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

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

Usage — Python (clean wrapper)

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

# 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

# list built-in schemes
prismalign schemes

Schemes

Python API

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)

Adapters (subprocess, prismalign.adapters)

All built on the shared CliAdapter base, which maps reads at three granularities:

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.

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; wfa2 is single-end for now (the subprocess sam/ strobealign/bowtie2/bwa_mem2 adapters are also SE).
• 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.