feat: paired-end alignment with adapter read-through detection

CHANGELOG.md

@@ -1,5 +1,37 @@
 # Changelog
 
+## Unreleased
+
+### Added
+
+- **Paired-end alignment** (#1). `-1 R1 -2 R2`, `--interleaved`, or a BAM of
+  0x1/0x40/0x80 mates align read pairs as pairs: the reference maximising the
+  two mates' summed score wins (both mates at or above `--min-score`), under an
+  orientation rule and `--max-insert` (default 500). A pair no reference
+  satisfies is written without 0x2, each mate at its own best reference on
+  either strand. Standard paired SAM: flags, `RNEXT`/`PNEXT`/`TLEN`, `MC`,
+  per-mate `NM`/`MD`/`AS`/`XS`/`XA`, MAPQ 60 only for a unique best reference,
+  mates adjacent in input order, `--sort coordinate` and `--secondary` supported.
+  See [docs/align.md](docs/align.md#paired-end).
+- 3' adapter read-through is detected before alignment (TruSeq Read 1/2
+  prefixes; `--adapter-r1`, `--adapter-r2`, `--min-adapter-overlap`,
+  `--no-adapter-trim`), the adapter is kept as soft clip with an `ad:i` tag, and
+  the rest of each mate is aligned end to end so terminal mismatches stay in
+  `MD`/`NM`. Paired runs therefore default to `--mode semiglobal` (an explicit
+  `--mode local` is honoured and warned about).
+- Library: `eschalign::adapter`, `Aligner::map_pairs`/`map_pairs_scored`,
+  `PairOptions`, `PairMapping`, `PairHit`.
+- Tests: `align_paired.rs` over a synthesized pair fixture in
+  `tests/fixtures/pairs/`; `align_unpaired.rs` pins the unpaired output to a
+  golden record-stream digest.
+
+### Changed
+
+- A BAM whose first record is flagged paired (0x1 with 0x40/0x80) is now a
+  paired run; it used to be read as unpaired reads with the flags discarded.
+  Every other input, and every existing flag, behaves as before (output
+  byte-identical, `unpaired_output_unchanged`).
+
 ## 0.1.0
 
 First release of eschalign as a standalone tool.

CLAUDE.md

@@ -183,13 +183,23 @@ would treat it as on. `clippy.toml`'s `disallowed-methods` enforces this (and
 
 - `crates/eschalign/` — alphabet / `Scoring` / `Mode` / scalar oracle / `simd` (two kernel families) / `sam` (CIGAR, `MD`/`NM`) / `panel` / `mapper` (`Aligner`) / `cuda` (`gpu`) / `fasta` / `error`. See its `CLAUDE.md` for the map, tests and traps.
 - `crates/eschalign-cli/src/` — `main.rs` (clap, tracing, rayon pool, allocator), `align.rs` (the pipeline), `device.rs` (`--device`), `env.rs`, `threads.rs`, `progress.rs`, `style.rs`, `uuid_util.rs`.
-- `crates/eschalign-cli/tests/` — `align_e2e.rs`, `align_sort.rs` over `tests/fixtures/` (60 real tRNA reads as a uBAM-shaped BAM plus the 47-reference panel and an ambiguity-code variant; **no pod5** — these tests never read signal).
+- `crates/eschalign-cli/tests/` — `align_e2e.rs`, `align_sort.rs`, `align_unpaired.rs` over `tests/fixtures/` (60 real tRNA reads as a uBAM-shaped BAM plus the 47-reference panel and an ambiguity-code variant; **no pod5** — these tests never read signal), and `align_paired.rs` over `tests/fixtures/pairs/` (a synthesized 4-reference panel with two paralogs and one hand-built read pair per behaviour; `gen_pair_fixture.py` regenerates it). `align_unpaired.rs` pins the unpaired output to a golden digest of the uncompressed record stream (`fixtures/unpaired_golden.tsv`, `ESCHALIGN_BLESS=1` rewrites it) — a change to it is a change to unpaired behaviour and must be intended.
 - `crates/eschalign/tests/fixtures/` — the parasail and calmd goldens and the scripts that regenerate them (`gen_align_golden.py` needs `parasail`; `gen_calmd_golden.py` needs samtools).
 
 ### The command
 
 `eschalign` — all-against-all alignment of reads to a small reference panel (tRNA) — `eschalign reads.ubam -r ref.fa -o out.bam`. Takes dorado's uBAM (or forward-strand aligned BAM, or FASTA/FASTQ, plain or gzip — sniffed from content), writes records in **input order** (`SO:unsorted`) — or, with `--sort coordinate` (#417), in exactly `samtools sort`'s order via a per-reference bucket sort in the writer thread (input order restored first, as the tie-break; past `--sort-memory` the largest buckets spill to one unlinked temp file, and each reference is ordered alone at the end, so there is no merge) — with **every input tag copied byte for byte** and `NM`/`MD`/`AS`/`XS`/`XA` owned by the aligner. Primary = lowest-indexed reference at the best score; the rest of the tie set goes to `XA` (bwa's `ref,±pos,CIGAR,NM;`) and, with `--secondary`, to 0x100 records carrying `SEQ`/`QUAL` as `*` and only the alignment tags + `RG`; MAPQ is 60 unique / 0 tied, nothing in between. Reads below `--min-score` are written unmapped, never dropped, and so are reads longer than `--max-read-len` (default 1000 nt, `0` = no limit; filtered in the CLI before any scoring, GPU included). Pipeline: reader thread (raw `bam::Record`s) → bounded channel → dispatcher cutting chunks of ≤ 64 reads / ≤ 16 kb and `spawn`ing each onto rayon (decode, score, trace the chunk's winners together, rewrite, **encode** — through a `bam::io::Writer` over a `Vec<u8>`, since `noodles_bam`'s record encoder is crate-private) → writer thread reordering numbered chunks into BGZF, returning one permit per chunk written (two batches in flight at most). The decode and encode live in the workers because a uBAM record is ~5 kB of `mv`/`MM`/`ML` for a ~130 nt read. **There is no barrier between batches, on evidence:** real tRNA samples carry reads up to ~400 kb, a chunk holding one costs seconds against milliseconds, and a `par_iter` per batch left 16 of 32 threads idle behind those stragglers (30.0 → 18.8 s on a 490 k-read sample, output identical); the per-batch pattern must not come back. `-v` logs the workers' busy time against wall × threads. **`--device` places `Stage::Align`** (#401): on the GPU, a dedicated thread sits between the reader and the dispatcher, takes whole batches off the reader's bounded channel, scores each in one `GpuScorer::score_batch` launch (forward reads, then reverse complements under `--strand both`), and hands the batch on with its `ScoreMatrix` behind an `Arc`; each chunk carries the `Arc` plus its offset, and its worker calls `Aligner::map_reads_scored` instead of `map_reads` — tie selection, pairs-kernel tracebacks, `MD`/`NM` and encoding are the CPU path's own code, so the output is byte-identical (`align_e2e::device_gpu_output_is_byte_identical`; also `samtools view | md5sum`-equal on a 490 k-read sample). A read the kernel does not take comes back unscored and its worker scores it on the CPU. The channel after the GPU thread holds two batches, so memory stays O(`--batch-size`). The GPU thread packs batch *k + 1* while batch *k*'s kernel runs (`GpuScorer::{prepare, submit, finish}`), so the device sits at ~100% whenever it has work; what leaves it idle is the CPU side stalling, and the very long reads (50–400 kb) are what stalls it (dropping the 52 reads over 20 kb takes M1 from 14.3 to 7.5 s): the GPU leaves them to the CPU, where the phase-1 score kernel is, as far as measured, **memory-bound on a read that long** (its per-row `H`/`E` state is 50 MB for a 395 kb read, streamed once per reference column: 4.9 s to score that read on one rna core, 2e9 cells/s), and the chunk holding it blocks the ordered writer while the permit window fills. Splitting such a read's scoring over the lane groups on six threads did not help *with that kernel* (six 50 MB streams at once), so #415 fixed the layout first: reads of ≥ 1,024 nt (`simd::TRANSPOSED_MIN_READ_LEN`) take a reference-major score kernel whose state is one panel row — 2.3× on the 395 kb read on one core — and a read of ≥ 16 kb (`SPLIT_MIN_BASES`) has its lane groups scored across the pool by `score_across_pool`, its row fed through the same `prescored` hook as a GPU row (0.29 s wall for that read; `ESCHALIGN_ROW_MAJOR_ONLY=1` is the A/B lever, output byte-identical). `--scoring` defaults to gpu-tRNA-mapper's `2,-1,-10,-1`; bwa `-x ont2d` is `1,-1,-2,-1`, and even then 27 of the 60 classify-fixture reads get a CIGAR different from bwa's (`tests/align_e2e.rs` prints the count; charging calls match the golden bit for bit on every read whose alignment is unchanged).
 
+### Paired-end (#1)
+
+`-1/-2`, `--interleaved`, or a BAM whose first primary record has 0x1 with 0x40/0x80 (`bam_is_paired` opens it once to peek) selects `Input::Paired` and a `PairCfg`; **nothing else changes for unpaired input** (`align_unpaired::unpaired_output_unchanged`). Where the pieces live:
+
+- **Library** (`eschalign`, still std only): `adapter.rs` (`find`, `detect_pair`: leftmost prefix match, `overlap/10` mismatches, the two mates' insert ends must agree or the weaker detection is dropped) and `mapper/pair.rs` (`Aligner::map_pairs_scored`, `PairOptions`, `PairMapping::{Proper, Separate}`). A *candidate* is `(reference, template strand)`; both mates must clear `min_score`; candidates are traced **tier by tier in descending summed score, every open pair's next tier through one `align_pairs` call** (so the traceback lanes still fill), and the first tier with a candidate that passes the constraints on the *traced* coordinates (no outie, insert ≤ `max_insert`, both mates align a base) is the tie set. Constraints need coordinates, so a score alone can never decide a pair. No qualifying candidate → each mate goes through `map_reads_scored` with `both_strands: true`.
+- **CLI** (`align.rs`): `PairSource` (`Bam` / `Interleaved` / `Files`) → `read_pair_batches` (whole pairs per batch, mate 1 then mate 2 adjacent, names stripped of `/1` `/2`, orphan/mismatch = error naming the read) → `cut_chunks(batch, paired)` (**a cut never falls between mates — the pair-in-one-chunk invariant; do not add a barrier either**) → `process_pair_chunk` → `build_pair_records` (flags, `RNEXT`/`PNEXT`/`TLEN`, `MC`, `ad`, per-mate `XA` entry *k* = tied pair *k*, secondaries per pair). The adapter bases are *not* aligned: the aligned mate is `seq[..len - ad]`, and `clipped_span` shifts a reverse hit's coordinates by `ad` so the soft clip lands at the right end.
+- **GPU**: the same thread, fed flattened reads (`2i` = mate 1, `2i+1` = mate 2) with `both_strands = true` always (the fallback needs every orientation). It must trim exactly what the workers trim — both call `pair_adapters`; if they ever diverge the rows are for different queries and nothing checks it but `pair_gpu_matches_cpu`.
+- **Paired default is `--mode semiglobal`**, because `local` clips the terminal mismatches the mode exists to keep; semi-global also gives every mate some positive score, which is why paired runs need `--min-score` (the tests all pass one).
+
+
 ## Testing notes
 
 - `cargo nextest run --workspace` is ~100 tests; the heaviest are the scalar-oracle equality tests (opt-level 2 in dev). `align_sort::sort_coordinate_matches_samtools_sort` uses `samtools` from PATH or `$SAMTOOLS` (in the pixi env; it says it skipped otherwise).

README.md

@@ -60,10 +60,32 @@ the content). Common options:
 | `--read-ids <FILE>` | align only the named reads |
 | `-t, --threads <N>`, `--batch-size <N>` | parallelism and memory |
 | `--device <auto\|cpu\|gpu>` | where panel scoring runs (see below) |
+| `-1`, `-2`, `--interleaved` | paired-end input — see [Paired-end](#paired-end) |
 | `-v`, `-vv`, `-q` | log verbosity; `RUST_LOG` wins |
 
 Full reference, output semantics and examples: [docs/align.md](docs/align.md).
 
+## Paired-end
+
+Illumina read pairs are aligned **as pairs**: the reference that maximises the
+two mates' summed score wins, so a mate that is ambiguous alone (paralogs,
+isodecoders) is placed by its partner. Input is two FASTQs, an interleaved
+FASTQ, or an unaligned BAM with 0x1/0x40/0x80 flags:
+
+```bash
+eschalign -1 R1.fq.gz -2 R2.fq.gz -r panel.fa -o aligned.bam --min-score 40
+eschalign pairs.fq.gz --interleaved -r panel.fa -o aligned.bam --min-score 40
+eschalign pairs.ubam -r panel.fa -o aligned.bam --min-score 40
+```
+
+3' adapter read-through is detected first (TruSeq prefixes by default) and kept
+as soft clip with an `ad:i` tag; the rest of each mate is aligned end to end, so
+a mismatch at a read end stays an `MD`/`NM` mismatch (paired runs default to
+`--mode semiglobal`). Pairs over `--max-insert` (default 500) or facing the wrong
+way are written without 0x2, each mate at its own best reference. Paired runs
+need a `--min-score`. See [docs/align.md](docs/align.md#paired-end) for the
+rules, flags and every decision.
+
 ## GPU
 
 A `--features gpu` build can score on a CUDA GPU (one thread per