diff --git a/CHANGELOG.md b/CHANGELOG.md index f9df0dd..6434479 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/CLAUDE.md b/CLAUDE.md index 385d94e..5743af3 100644 --- a/CLAUDE.md +++ b/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`, 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). diff --git a/README.md b/README.md index 45daa9e..5ef3981 100644 --- a/README.md +++ b/README.md @@ -60,10 +60,32 @@ the content). Common options: | `--read-ids ` | align only the named reads | | `-t, --threads `, `--batch-size ` | parallelism and memory | | `--device ` | 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 diff --git a/crates/eschalign-cli/src/align.rs b/crates/eschalign-cli/src/align.rs index 4a2f862..8e0f21b 100644 --- a/crates/eschalign-cli/src/align.rs +++ b/crates/eschalign-cli/src/align.rs @@ -120,8 +120,8 @@ use sam::header::record::value::map::{self, Program, ReferenceSequence, program: use std::sync::Arc; use eschalign::{ - Aligner, Backend, CigarOp, Hit, MapOptions, Mode, Panel, ReadMapping, ScoreMatrix, Scoring, - sam as align_sam, + Aligner, Backend, CigarOp, Hit, MapOptions, Mode, PairMapping, PairOptions, Panel, ReadMapping, + ScoreMatrix, Scoring, sam as align_sam, }; use rayon::prelude::*; @@ -133,8 +133,50 @@ pub struct AlignArgs { /// Reads: an unaligned BAM (dorado's uBAM), or FASTA/FASTQ, plain or gzip. /// Forward-strand aligned BAM records are accepted too; their alignment is /// discarded and they are aligned afresh - #[arg(value_name = "READS")] - pub reads: PathBuf, + #[arg( + value_name = "READS", + required_unless_present_all = ["r1", "r2"], + conflicts_with_all = ["r1", "r2"] + )] + pub reads: Option, + + /// Mate 1 FASTQ/FASTA (plain or gzip) of a read pair; needs `-2`. Mates + /// are matched by name, `/1` and `/2` stripped. Switches on paired-end + /// alignment + #[arg(short = '1', long = "r1", value_name = "FASTQ", requires = "r2")] + pub r1: Option, + + /// Mate 2 FASTQ/FASTA of a read pair; needs `-1` + #[arg(short = '2', long = "r2", value_name = "FASTQ", requires = "r1")] + pub r2: Option, + + /// READS is an interleaved FASTQ/FASTA: mate 1, mate 2, mate 1, … (a BAM's + /// pairing is read from its 0x1/0x40/0x80 flags instead) + #[arg(long, requires = "reads", conflicts_with_all = ["r1", "r2"])] + pub interleaved: bool, + + /// Paired input only: the longest insert (leftmost to rightmost aligned + /// reference base of either mate) a proper pair may have [default: 500] + #[arg(long, value_name = "N")] + pub max_insert: Option, + + /// Paired input only: mate 1's 3' adapter (its prefix is enough) + /// [default: TruSeq Read 1, AGATCGGAAGAGCACACGTCTGAACTCCAGTCA] + #[arg(long, value_name = "SEQ")] + pub adapter_r1: Option, + + /// Paired input only: mate 2's 3' adapter + /// [default: TruSeq Read 2, AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT] + #[arg(long, value_name = "SEQ")] + pub adapter_r2: Option, + + /// Paired input only: the shortest 3' match counted as adapter [default: 5] + #[arg(long, value_name = "N")] + pub min_adapter_overlap: Option, + + /// Paired input only: do not look for 3' adapter read-through (no `ad` tag) + #[arg(long)] + pub no_adapter_trim: bool, /// Reference panel FASTA, plain or gzip #[arg(short, long, value_name = "FASTA")] @@ -166,9 +208,10 @@ pub struct AlignArgs { pub tmp_dir: Option, /// `local` (Smith-Waterman, overhangs soft-clipped) or `semiglobal` - /// (overlap: leading/trailing gaps of either sequence free) - #[arg(long, default_value = "local", value_parser = parse_mode, value_name = "MODE")] - pub mode: Mode, + /// (overlap: leading/trailing gaps of either sequence free) [default: + /// `local`; `semiglobal` for paired input, which keeps terminal mismatches] + #[arg(long, value_parser = parse_mode, value_name = "MODE")] + pub mode: Option, /// MATCH,MISMATCH,GAP_OPEN,GAP_EXTEND. A gap of length k scores /// GAP_OPEN + (k-1)*GAP_EXTEND. bwa's `-A1 -B1 -O1 -E1` is `1,-1,-2,-1` @@ -612,6 +655,323 @@ enum Source { Fastx(FastxReader), } +/// What the reader thread will read: single reads, or read pairs. +enum Input { + Single(Source), + Paired(PairSource), +} + +// --------------------------------------------------------------------------- +// Paired input +// --------------------------------------------------------------------------- + +/// Where read pairs come from. Whichever it is, the reader hands over pairs +/// as two adjacent records, mate 1 then mate 2. +enum PairSource { + /// A BAM whose records carry 0x1 with 0x40 / 0x80, mates adjacent. + Bam(bam::io::Reader>), + /// One FASTQ/FASTA, mate 1 then mate 2 for each pair. + Interleaved(FastxReader), + /// Two FASTQ/FASTA files, record `k` of each being a pair. + Files(Box, Box, [String; 2]), +} + +/// A read name without a trailing `/1` or `/2` (the name both mates share). +fn mate_stem(name: &[u8]) -> &[u8] { + match name { + [stem @ .., b'/', b'1' | b'2'] => stem, + _ => name, + } +} + +fn show(name: &[u8]) -> String { + String::from_utf8_lossy(name).into_owned() +} + +/// Whether a BAM holds read pairs: its first primary record is flagged +/// paired (0x1) and as mate 1 or mate 2 (0x40 / 0x80). +fn bam_is_paired(path: &Path) -> anyhow::Result { + let f = File::open(path).with_context(|| format!("opening {}", path.display()))?; + let mut r = bam::io::Reader::new(f); + r.read_header() + .with_context(|| format!("reading BAM header of {}", path.display()))?; + let mut rec = bam::Record::default(); + while r.read_record(&mut rec)? != 0 { + let flags = rec.flags(); + if flags.is_secondary() || flags.is_supplementary() { + continue; + } + return Ok(flags.is_segmented() && (flags.is_first_segment() || flags.is_last_segment())); + } + Ok(false) +} + +/// The next primary record of a BAM, with the same refusals as the unpaired +/// reader (secondary/supplementary skipped and counted, reverse-strand +/// alignments refused). +fn next_primary( + r: &mut bam::io::Reader>, + counts: &mut ReadCounts, +) -> anyhow::Result> { + loop { + let mut rec = bam::Record::default(); + if r.read_record(&mut rec)? == 0 { + return Ok(None); + } + let flags = rec.flags(); + if flags.is_secondary() || flags.is_supplementary() { + counts.not_primary += 1; + continue; + } + if !flags.is_unmapped() && flags.is_reverse_complemented() { + bail!( + "input record {} is aligned to the reverse strand; `eschalign` takes \ + unaligned reads (or forward-strand alignments, whose SEQ is the read as \ + sequenced)", + rec.name() + .map(|n| n.to_string()) + .unwrap_or_else(|| "*".into()) + ); + } + return Ok(Some(rec)); + } +} + +impl PairSource { + /// The next pair as `(mate 1, mate 2, shared name)`, or `None` at the end. + /// An orphan mate or a name mismatch is an error naming the read. + fn next_pair( + &mut self, + counts: &mut ReadCounts, + ) -> anyhow::Result)>> { + match self { + PairSource::Bam(r) => { + let Some(a) = next_primary(r, counts)? else { + return Ok(None); + }; + let name_a = a.name().map(|n| n.to_vec()).unwrap_or_default(); + let fa = a.flags(); + if !fa.is_segmented() || fa.is_first_segment() == fa.is_last_segment() { + bail!( + "read {}: not a mate (flags {}): paired input needs 0x1 with exactly \ + one of 0x40 / 0x80 on every record", + show(&name_a), + u16::from(fa) + ); + } + let Some(b) = next_primary(r, counts)? else { + bail!("read {} has no mate: the BAM ends after it", show(&name_a)); + }; + let name_b = b.name().map(|n| n.to_vec()).unwrap_or_default(); + let fb = b.flags(); + if mate_stem(&name_a) != mate_stem(&name_b) { + bail!( + "read {} has no mate: the next record is {} (mates must be adjacent, \ + name-grouped)", + show(&name_a), + show(&name_b) + ); + } + if !fb.is_segmented() + || fb.is_first_segment() == fb.is_last_segment() + || fb.is_first_segment() == fa.is_first_segment() + { + bail!( + "read {}: the two records are not mate 1 and mate 2 (flags {} and {})", + show(mate_stem(&name_a)), + u16::from(fa), + u16::from(fb) + ); + } + let stem = mate_stem(&name_a).to_vec(); + let (r1, r2) = if fa.is_first_segment() { + (a, b) + } else { + (b, a) + }; + Ok(Some((InRecord::Bam(r1), InRecord::Bam(r2), stem))) + } + PairSource::Interleaved(r) => { + let Some((na, sa, qa)) = r.next_record()? else { + return Ok(None); + }; + let Some((nb, sb, qb)) = r.next_record()? else { + bail!("read {} has no mate: the file ends after it", show(&na)); + }; + if mate_stem(&na) != mate_stem(&nb) { + bail!( + "read {} has no mate: the next record is {} (an interleaved file \ + alternates mate 1, mate 2)", + show(&na), + show(&nb) + ); + } + let stem = mate_stem(&na).to_vec(); + Ok(Some(( + InRecord::Buf(fastx_to_record(stem.clone(), sa, qa)), + InRecord::Buf(fastx_to_record(stem.clone(), sb, qb)), + stem, + ))) + } + PairSource::Files(f1, f2, names) => match (f1.next_record()?, f2.next_record()?) { + (None, None) => Ok(None), + (Some((n, ..)), None) => bail!( + "read {} has no mate: {} ends before {}", + show(&n), + names[1], + names[0] + ), + (None, Some((n, ..))) => bail!( + "read {} has no mate: {} ends before {}", + show(&n), + names[0], + names[1] + ), + (Some((na, sa, qa)), Some((nb, sb, qb))) => { + if mate_stem(&na) != mate_stem(&nb) { + bail!( + "mate names differ: {} in {} but {} in {} at the same record", + show(&na), + names[0], + show(&nb), + names[1] + ); + } + let stem = mate_stem(&na).to_vec(); + Ok(Some(( + InRecord::Buf(fastx_to_record(stem.clone(), sa, qa)), + InRecord::Buf(fastx_to_record(stem.clone(), sb, qb)), + stem, + ))) + } + }, + } + } +} + +/// The reader thread's body for pairs: batches of whole pairs (mate 1, mate 2 +/// adjacent) until the end. `--read-ids` selects by pair name. +fn read_pair_batches( + mut source: PairSource, + wanted: Option>>, + batch_size: usize, + tx: SyncSender>, +) -> anyhow::Result { + let mut counts = ReadCounts::default(); + let mut batch = Vec::with_capacity(batch_size + 1); + while let Some((a, b, name)) = source.next_pair(&mut counts)? { + if let Some(w) = &wanted + && !w.contains(&canonical_read_name(&name)) + { + counts.filtered += 2; + continue; + } + batch.push(a); + batch.push(b); + if batch.len() >= batch_size { + let full = std::mem::replace(&mut batch, Vec::with_capacity(batch_size + 1)); + if tx.send(full).is_err() { + return Ok(counts); + } + } + } + if !batch.is_empty() { + let _ = tx.send(batch); + } + Ok(counts) +} + +/// Paired-run settings, shared by the workers and the GPU thread. +#[derive(Clone)] +struct PairCfg { + max_insert: usize, + /// `None` under `--no-adapter-trim`. + adapters: Option, +} + +#[derive(Clone)] +struct AdapterSet { + r1: Vec, + r2: Vec, + min_overlap: usize, +} + +impl AdapterSet { + fn config(&self) -> eschalign::adapter::AdapterConfig<'_> { + eschalign::adapter::AdapterConfig { + r1: &self.r1, + r2: &self.r2, + min_overlap: self.min_overlap, + } + } +} + +/// Adapter bases at the 3' end of each mate of a pair. A mate over +/// `--max-read-len` is not aligned and is not searched (0); the other mate is +/// searched alone, since the two-mate agreement check needs both. One function +/// for the workers and the GPU thread, so both trim the same bases. +fn pair_adapters(l1: &[u8], l2: &[u8], cfg: &PairCfg, max_read_len: usize) -> [usize; 2] { + let Some(set) = &cfg.adapters else { + return [0, 0]; + }; + let (o1, o2) = ( + over_limit(max_read_len, l1.len()), + over_limit(max_read_len, l2.len()), + ); + let ac = set.config(); + match (o1, o2) { + (false, false) => eschalign::adapter::detect_pair(l1, l2, &ac), + (false, true) => [eschalign::adapter::find(l1, ac.r1, ac.min_overlap), 0], + (true, false) => [0, eschalign::adapter::find(l2, ac.r2, ac.min_overlap)], + (true, true) => [0, 0], + } +} + +/// The paired-run settings from the command line. +fn pair_cfg_from(args: &AlignArgs) -> anyhow::Result { + use eschalign::adapter::{DEFAULT_MIN_OVERLAP, TRUSEQ_R1, TRUSEQ_R2}; + let adapters = if args.no_adapter_trim { + if args.adapter_r1.is_some() + || args.adapter_r2.is_some() + || args.min_adapter_overlap.is_some() + { + bail!("--no-adapter-trim contradicts --adapter-r1/--adapter-r2/--min-adapter-overlap"); + } + None + } else { + let r1 = match &args.adapter_r1 { + Some(a) => parse_adapter("--adapter-r1", a)?, + None => TRUSEQ_R1.to_vec(), + }; + let r2 = match &args.adapter_r2 { + Some(a) => parse_adapter("--adapter-r2", a)?, + None => TRUSEQ_R2.to_vec(), + }; + let min_overlap = args.min_adapter_overlap.unwrap_or(DEFAULT_MIN_OVERLAP); + if min_overlap == 0 { + bail!("--min-adapter-overlap must be at least 1"); + } + Some(AdapterSet { + r1, + r2, + min_overlap, + }) + }; + Ok(PairCfg { + max_insert: args.max_insert.unwrap_or(500), + adapters, + }) +} + +/// One adapter sequence from the command line: A/C/G/T/U/N, any case. +fn parse_adapter(flag: &str, s: &str) -> anyhow::Result> { + let seq: Vec = s.trim().bytes().map(|b| b.to_ascii_uppercase()).collect(); + if seq.is_empty() || !seq.iter().all(|b| b"ACGTUN".contains(b)) { + bail!("{flag}: {s:?} is not a nucleotide sequence (A, C, G, T/U, N)"); + } + Ok(seq) +} + // --------------------------------------------------------------------------- // Output records // --------------------------------------------------------------------------- @@ -634,10 +994,26 @@ fn int_value(v: i64) -> Value { } fn sam_cigar(hit: &Hit, read_len: usize) -> sam::alignment::record_buf::Cigar { + sam_cigar_ad(hit, read_len, 0) +} + +/// Where a hit's aligned part sits in the full SEQ as written: the alignment +/// is on the mate *without* its `ad` adapter bases, which are the mate's 3' +/// end — the end of SEQ for a forward hit, the start for a reverse one. +fn clipped_span(hit: &Hit, ad: usize) -> (usize, usize) { + let a = &hit.alignment; + let shift = if hit.reverse { ad } else { 0 }; + (a.query_start + shift, a.query_end + shift) +} + +/// [`sam_cigar`] for a read whose last `ad` 3' bases were not aligned and are +/// soft-clipped (`ad = 0`: the plain unpaired CIGAR). +fn sam_cigar_ad(hit: &Hit, read_len: usize, ad: usize) -> sam::alignment::record_buf::Cigar { let a = &hit.alignment; + let (qs, qe) = clipped_span(hit, ad); let mut ops = Vec::with_capacity(a.ops.len() + 2); - if a.query_start > 0 { - ops.push(Op::new(Kind::SoftClip, a.query_start)); + if qs > 0 { + ops.push(Op::new(Kind::SoftClip, qs)); } for &(op, len) in &a.ops { let kind = match op { @@ -647,21 +1023,26 @@ fn sam_cigar(hit: &Hit, read_len: usize) -> sam::alignment::record_buf::Cigar { }; ops.push(Op::new(kind, len as usize)); } - if read_len > a.query_end { - ops.push(Op::new(Kind::SoftClip, read_len - a.query_end)); + if read_len > qe { + ops.push(Op::new(Kind::SoftClip, read_len - qe)); } ops.into() } /// bwa's `XA:Z:` entry: `ref,±pos,CIGAR,NM;`. fn xa_entry(panel: &Panel, hit: &Hit, read_len: usize) -> String { + xa_entry_ad(panel, hit, read_len, 0) +} + +fn xa_entry_ad(panel: &Panel, hit: &Hit, read_len: usize, ad: usize) -> String { let a = &hit.alignment; + let (qs, qe) = clipped_span(hit, ad); format!( "{},{}{},{},{};", panel.get(hit.reference).name, if hit.reverse { '-' } else { '+' }, a.ref_start + 1, - align_sam::cigar_string(&a.ops, read_len, a.query_start, a.query_end), + align_sam::cigar_string(&a.ops, read_len, qs, qe), hit.nm ) } @@ -781,6 +1162,230 @@ fn build_records( ) } +/// Tags a paired run owns besides [`OWNED_TAGS`]: the mate's CIGAR and the +/// adapter length. +const MC_TAG: Tag = Tag::new(b'M', b'C'); +const AD_TAG: Tag = Tag::new(b'a', b'd'); + +/// How one mate was placed, however its pair resolved. +struct MatePlacement<'a> { + /// `None`: the mate is unmapped. + primary: Option<&'a Hit>, + /// The mate's hit in every other tie-set member (`XA`, `--secondary`). + others: Vec<&'a Hit>, + /// Beside each of `others`, the partner's hit it belongs to: in a proper + /// pair the other mate at the same tie-set member, otherwise the + /// partner's primary. + partners: Vec>, + /// `AS` of the primary. + score: i32, + xs: Option, + /// MAPQ 60 (a unique best reference) rather than 0. + unique: bool, +} + +/// `TLEN` of `hit` against its partner: the span from the leftmost aligned +/// base of either to the rightmost, positive for the leftmost mate (mate 1 +/// when both start together) and negative for the other. +fn template_length(hit: &Hit, partner: &Hit, first: bool) -> i32 { + let (a, b) = (&hit.alignment, &partner.alignment); + let len = a.ref_end.max(b.ref_end) - a.ref_start.min(b.ref_start); + let left = a.ref_start < b.ref_start || (a.ref_start == b.ref_start && first); + let len = i32::try_from(len).unwrap_or(i32::MAX); + if left { len } else { -len } +} + +/// Rewrite a read pair's two input records into its output records: the two +/// primaries (mate 1, mate 2) and, with `--secondary`, one secondary per mate +/// per further tie-set member. `ad` is `Some` when adapters were searched for +/// (it carries each mate's `ad:i`). +fn build_pair_records( + mut recs: [RecordBuf; 2], + mates: &[MatePlacement<'_>; 2], + proper: bool, + ad: Option<[usize; 2]>, + panel: &Panel, + secondary: bool, +) -> ([RecordBuf; 2], Vec, [Outcome; 2]) { + let full = [recs[0].sequence().len(), recs[1].sequence().len()]; + let adn = ad.unwrap_or([0, 0]); + let seg = [Flags::FIRST_SEGMENT, Flags::LAST_SEGMENT]; + let mut flags = [Flags::empty(); 2]; + let mut kept = [Flags::empty(); 2]; + + // Each mate on its own: strand, position, CIGAR, alignment tags. + for m in 0..2 { + let rec = &mut recs[m]; + kept[m] = rec.flags() & (Flags::QC_FAIL | Flags::DUPLICATE); + for t in OWNED_TAGS.iter().chain([&MC_TAG, &AD_TAG]) { + rec.data_mut().remove(t); + } + *rec.mate_reference_sequence_id_mut() = None; + *rec.mate_alignment_start_mut() = None; + *rec.template_length_mut() = 0; + flags[m] = kept[m] | Flags::SEGMENTED | seg[m]; + let place = &mates[m]; + match place.primary { + None => { + flags[m] |= Flags::UNMAPPED; + *rec.reference_sequence_id_mut() = None; + *rec.alignment_start_mut() = None; + *rec.mapping_quality_mut() = MappingQuality::new(0); + *rec.cigar_mut() = Default::default(); + } + Some(h) => { + if h.reverse { + flags[m] |= Flags::REVERSE_COMPLEMENTED; + let rc = eschalign::alphabet::reverse_complement(rec.sequence().as_ref()); + *rec.sequence_mut() = rc.into(); + let mut q: Vec = rec.quality_scores().as_ref().to_vec(); + q.reverse(); + *rec.quality_scores_mut() = q.into(); + } + *rec.reference_sequence_id_mut() = Some(h.reference); + *rec.alignment_start_mut() = Some(position(h.alignment.ref_start)); + *rec.mapping_quality_mut() = MappingQuality::new(if place.unique { 60 } else { 0 }); + *rec.cigar_mut() = sam_cigar_ad(h, full[m], adn[m]); + let data = rec.data_mut(); + data.insert(Tag::EDIT_DISTANCE, int_value(h.nm as i64)); + data.insert( + Tag::MISMATCHED_POSITIONS, + Value::String(h.md.clone().into()), + ); + data.insert(Tag::ALIGNMENT_SCORE, int_value(place.score as i64)); + if let Some(xs) = place.xs { + data.insert(Tag::new(b'X', b'S'), int_value(xs as i64)); + } + if !place.others.is_empty() { + let xa: String = place + .others + .iter() + .map(|o| xa_entry_ad(panel, o, full[m], adn[m])) + .collect(); + data.insert(Tag::new(b'X', b'A'), Value::String(xa.into())); + } + } + } + if let Some(a) = ad { + rec.data_mut().insert(AD_TAG, int_value(a[m] as i64)); + } + } + + // Each mate against its partner: the mate fields. + let primaries = [mates[0].primary, mates[1].primary]; + for m in 0..2 { + let o = 1 - m; + let rec = &mut recs[m]; + match (primaries[m], primaries[o]) { + (Some(h), Some(oh)) => { + if oh.reverse { + flags[m] |= Flags::MATE_REVERSE_COMPLEMENTED; + } + if proper { + flags[m] |= Flags::PROPERLY_SEGMENTED; + } + *rec.mate_reference_sequence_id_mut() = Some(oh.reference); + *rec.mate_alignment_start_mut() = Some(position(oh.alignment.ref_start)); + let mc = align_sam::cigar_string( + &oh.alignment.ops, + full[o], + clipped_span(oh, adn[o]).0, + clipped_span(oh, adn[o]).1, + ); + rec.data_mut().insert(MC_TAG, Value::String(mc.into())); + if h.reference == oh.reference { + *rec.template_length_mut() = template_length(h, oh, m == 0); + } + } + (Some(h), None) => { + flags[m] |= Flags::MATE_UNMAPPED; + *rec.mate_reference_sequence_id_mut() = Some(h.reference); + *rec.mate_alignment_start_mut() = Some(position(h.alignment.ref_start)); + } + (None, Some(oh)) => { + // An unmapped mate sits at its partner's position, as samtools + // writes it. + if oh.reverse { + flags[m] |= Flags::MATE_REVERSE_COMPLEMENTED; + } + *rec.reference_sequence_id_mut() = Some(oh.reference); + *rec.alignment_start_mut() = Some(position(oh.alignment.ref_start)); + *rec.mate_reference_sequence_id_mut() = Some(oh.reference); + *rec.mate_alignment_start_mut() = Some(position(oh.alignment.ref_start)); + } + (None, None) => flags[m] |= Flags::MATE_UNMAPPED, + } + *rec.flags_mut() = flags[m]; + } + + // --secondary: each further tie-set member as a pair of 0x100 records, + // mate 1's then mate 2's, in tie order. + let mut secs = Vec::new(); + if secondary { + let n_ties = mates[0].others.len().max(mates[1].others.len()); + for k in 0..n_ties { + for m in 0..2 { + let Some(&h) = mates[m].others.get(k) else { + continue; + }; + let partner = mates[m].partners[k]; + let mut f = Flags::SECONDARY | Flags::SEGMENTED | seg[m] | kept[m]; + if h.reverse { + f |= Flags::REVERSE_COMPLEMENTED; + } + let mut d = sam::alignment::record_buf::Data::default(); + d.insert(Tag::EDIT_DISTANCE, int_value(h.nm as i64)); + d.insert( + Tag::MISMATCHED_POSITIONS, + Value::String(h.md.clone().into()), + ); + d.insert(Tag::ALIGNMENT_SCORE, int_value(h.alignment.score as i64)); + if let Some(rg) = recs[m].data().get(&Tag::READ_GROUP) { + d.insert(Tag::READ_GROUP, rg.clone()); + } + let mut b = RecordBuf::builder() + .set_reference_sequence_id(h.reference) + .set_alignment_start(position(h.alignment.ref_start)) + .set_mapping_quality(MappingQuality::new(0).expect("0 is a valid MAPQ")) + .set_cigar(sam_cigar_ad(h, full[m], adn[m])); + match partner { + Some(ph) => { + if ph.reverse { + f |= Flags::MATE_REVERSE_COMPLEMENTED; + } + b = b + .set_mate_reference_sequence_id(ph.reference) + .set_mate_alignment_start(position(ph.alignment.ref_start)); + if proper { + f |= Flags::PROPERLY_SEGMENTED; + if ph.reference == h.reference { + b = b.set_template_length(template_length(h, ph, m == 0)); + } + } + } + None => { + f |= Flags::MATE_UNMAPPED; + b = b + .set_mate_reference_sequence_id(h.reference) + .set_mate_alignment_start(position(h.alignment.ref_start)); + } + } + b = b.set_flags(f).set_data(d); + if let Some(n) = recs[m].name() { + b = b.set_name(n.to_vec()); + } + secs.push(b.build()); + } + } + } + let outcome = |m: usize| match mates[m].primary { + None => Outcome::Unmapped, + Some(_) if mates[m].unique => Outcome::Unique, + Some(_) => Outcome::Tied, + }; + (recs, secs, [outcome(0), outcome(1)]) +} + /// Everything a worker needs, shared read-only across the pool. struct Job<'a> { aligner: &'a Aligner, @@ -797,6 +1402,11 @@ struct Job<'a> { max_read_len: usize, /// Reads written unmapped for `--max-read-len`. skipped_long: AtomicU64, + /// `Some` for a paired run: the chunks hold read pairs. + pair: Option, + /// Read pairs seen, and of those how many resolved as proper pairs. + pairs: AtomicU64, + proper_pairs: AtomicU64, } /// Whether `--max-read-len` keeps a read of `len` bases out of alignment. @@ -937,6 +1547,136 @@ fn process_chunk( .collect() } +/// A chunk of read pairs (mate 1, mate 2 adjacent), start to finish on a +/// worker: decode, find adapters, align each pair as a pair, rewrite, encode. +fn process_pair_chunk( + ctx: &Job<'_>, + cfg: &PairCfg, + chunk: Vec, + scores: ChunkScores, +) -> anyhow::Result { + let mut records = chunk + .into_iter() + .map(|rec| match rec { + InRecord::Bam(raw) => Ok(RecordBuf::try_from_alignment_record(ctx.in_header, &raw)?), + InRecord::Buf(r) => Ok(r), + }) + .collect::>>()?; + anyhow::ensure!( + records.len().is_multiple_of(2), + "a chunk of pairs has an odd record count" + ); + // Both mates share one name, without any /1 or /2. + for r in &mut records { + if let Some(n) = r.name() { + let stem = mate_stem(n.as_ref()); + if stem.len() != n.len() { + *r.name_mut() = Some(stem.to_vec().into()); + } + } + } + let n = records.len() / 2; + let max = ctx.max_read_len; + let seq = |k: usize| -> &[u8] { records[k].sequence().as_ref() }; + let ads: Vec<[usize; 2]> = (0..n) + .map(|i| pair_adapters(seq(2 * i), seq(2 * i + 1), cfg, max)) + .collect(); + // What is aligned: each mate less its 3' adapter; a mate over + // --max-read-len is not aligned at all (and the other is mapped alone). + let mut skipped = 0u64; + let trimmed: Vec<[&[u8]; 2]> = (0..n) + .map(|i| { + [0, 1].map(|m| { + let s = seq(2 * i + m); + if over_limit(max, s.len()) { + skipped += 1; + &s[..0] + } else { + &s[..s.len() - ads[i][m]] + } + }) + }) + .collect(); + if skipped > 0 { + ctx.skipped_long.fetch_add(skipped, Ordering::Relaxed); + } + let popts = PairOptions { + min_score: ctx.opts.min_score, + max_ties: ctx.opts.max_ties, + max_insert: cfg.max_insert, + both_strands: ctx.opts.both_strands, + }; + let mappings = ctx + .aligner + .map_pairs_scored(&trimmed, &popts, |i, m, reverse| { + scores + .as_ref() + .and_then(|(batch, first)| batch.row(first + 2 * i + m, reverse)) + }); + drop(trimmed); + ctx.pairs.fetch_add(n as u64, Ordering::Relaxed); + let mut out: Encoded = Vec::with_capacity(records.len()); + let mut recs = records.into_iter(); + for (mapping, ad) in mappings.iter().zip(&ads) { + let pair = [recs.next().expect("mate 1"), recs.next().expect("mate 2")]; + let (placements, proper) = match mapping { + PairMapping::Proper { + n_tied, + suboptimal, + hits, + .. + } => { + let place = |m: usize| MatePlacement { + primary: Some(&hits[0].mates[m]), + others: hits[1..].iter().map(|h| &h.mates[m]).collect(), + partners: hits[1..].iter().map(|h| Some(&h.mates[1 - m])).collect(), + score: hits[0].mates[m].alignment.score, + xs: suboptimal[m], + unique: *n_tied == 1, + }; + ctx.proper_pairs.fetch_add(1, Ordering::Relaxed); + ([place(0), place(1)], true) + } + PairMapping::Separate(ms) => { + let place = |m: usize| { + let me = &ms[m]; + let others: Vec<&Hit> = me.hits.iter().skip(1).collect(); + MatePlacement { + partners: vec![ms[1 - m].primary(); others.len()], + others, + primary: me.primary(), + score: me.best_score.unwrap_or(0), + xs: me.suboptimal, + unique: me.n_tied == 1, + } + }; + ([place(0), place(1)], false) + } + }; + let (prim, secs, outcomes) = build_pair_records( + pair, + &placements, + proper, + cfg.adapters.is_some().then_some(*ad), + ctx.aligner.panel(), + ctx.secondary, + ); + let encode = |recs: &[&RecordBuf]| -> anyhow::Result> { + let mut w = bam::io::Writer::from(Vec::with_capacity(8 << 10)); + for r in recs { + w.write_alignment_record(ctx.out_header, *r)?; + } + Ok(w.into_inner()) + }; + let [r1, r2] = &prim; + out.push((encode(&[r1])?, outcomes[0])); + let mut tail: Vec<&RecordBuf> = vec![r2]; + tail.extend(secs.iter()); + out.push((encode(&tail)?, outcomes[1])); + } + Ok(out) +} + /// One read strand's panel scores (codes in), computed one lane group per /// rayon task — the row `Aligner::score_all` would give, as `i16` for /// `map_reads_scored`. `None` outside the `i16` bound, where the aligner's @@ -1022,6 +1762,120 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { bail!("--batch-size must be at least 1"); } + // --- Input ------------------------------------------------------------- + let wanted = args.read_ids.as_deref().map(load_read_ids).transpose()?; + if let Some(w) = &wanted { + info!( + "aligning only the {} reads named in --read-ids", + style::count(w.len()) + ); + } + let threads = rayon::current_num_threads().max(1); + let open_bam = |path: &Path| -> anyhow::Result<(_, sam::Header)> { + let workers = std::num::NonZero::new(threads.div_ceil(4)).expect("at least one"); + let file = File::open(path)?; + let mut r = bam::io::Reader::from(bgzf::io::MultithreadedReader::with_worker_count( + workers, file, + )); + let h = r + .read_header() + .with_context(|| format!("reading BAM header of {}", path.display()))?; + Ok((r, h)) + }; + let open_fastx = |path: &Path| -> anyhow::Result { + match sniff(path)? { + InputFormat::Fasta => Ok(FastxReader::new(open_text(path)?, false)), + InputFormat::Fastq => Ok(FastxReader::new(open_text(path)?, true)), + InputFormat::Bam => bail!( + "{}: a BAM cannot be a mate file; pass a paired BAM as READS", + path.display() + ), + } + }; + let mut in_header = sam::Header::default(); + let input = if let (Some(r1), Some(r2)) = (&args.r1, &args.r2) { + let names = [r1.display().to_string(), r2.display().to_string()]; + Input::Paired(PairSource::Files( + Box::new(open_fastx(r1)?), + Box::new(open_fastx(r2)?), + names, + )) + } else { + let reads = args.reads.as_deref().expect("clap requires READS or -1/-2"); + match sniff(reads)? { + InputFormat::Bam => { + if args.interleaved { + bail!( + "--interleaved applies to FASTA/FASTQ; a BAM's pairing is read from \ + its 0x1/0x40/0x80 flags" + ); + } + let paired = bam_is_paired(reads)?; + let (r, h) = open_bam(reads)?; + in_header = h; + if paired { + Input::Paired(PairSource::Bam(r)) + } else { + Input::Single(Source::Bam(r)) + } + } + InputFormat::Fasta | InputFormat::Fastq if args.interleaved => { + Input::Paired(PairSource::Interleaved(open_fastx(reads)?)) + } + InputFormat::Fasta => { + Input::Single(Source::Fastx(FastxReader::new(open_text(reads)?, false))) + } + InputFormat::Fastq => { + Input::Single(Source::Fastx(FastxReader::new(open_text(reads)?, true))) + } + } + }; + let paired = matches!(input, Input::Paired(_)); + let pair_cfg = if paired { + Some(pair_cfg_from(&args)?) + } else { + let given = [ + ("--max-insert", args.max_insert.is_some()), + ("--adapter-r1", args.adapter_r1.is_some()), + ("--adapter-r2", args.adapter_r2.is_some()), + ("--min-adapter-overlap", args.min_adapter_overlap.is_some()), + ("--no-adapter-trim", args.no_adapter_trim), + ]; + let named: Vec<&str> = given.iter().filter(|g| g.1).map(|g| g.0).collect(); + if !named.is_empty() { + warn!( + "{} only apply to paired input (-1/-2, --interleaved, or a BAM of 0x40/0x80 \ + mates); this input is unpaired, so they are ignored", + named.join(", ") + ); + } + None + }; + // Paired input aligns end to end on the reference by default, so a + // mismatch at a read end stays a mismatch; `local` would clip it. + let mode = args.mode.unwrap_or(if paired { + Mode::SemiGlobal + } else { + Mode::Local + }); + if paired { + info!( + "paired-end: {} mode, max insert {}, adapters {}", + mode.name(), + pair_cfg.as_ref().map_or(0, |c| c.max_insert), + match pair_cfg.as_ref().and_then(|c| c.adapters.as_ref()) { + Some(a) => format!("searched (min overlap {})", a.min_overlap), + None => "not searched (--no-adapter-trim)".to_string(), + } + ); + if mode == Mode::Local { + warn!( + "--mode local soft-clips mismatches near a read end, which erases the \ + misincorporation and RT-stop signal paired mode is meant to keep" + ); + } + } + // --- Reference panel + kernel ------------------------------------------ let panel = load_reference(&args.reference)?; info!( @@ -1046,7 +1900,7 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { } let row_major_only = crate::env::flag_with_alias(ROW_MAJOR_ONLY_ENV, ROW_MAJOR_ONLY_ENV_DEPRECATED); - let mut aligner = Aligner::new(panel, args.scoring, args.mode)?; + let mut aligner = Aligner::new(panel, args.scoring, mode)?; if row_major_only { aligner = aligner.with_transposed_min_len(None); info!( @@ -1056,7 +1910,7 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { } info!( "{} mode, scoring {} (match,mismatch,gap_open,gap_extend), {} strand{}, kernel {} ({} lanes)", - args.mode.name(), + mode.name(), args.scoring, match args.strand { StrandArg::Forward => "forward", @@ -1081,36 +1935,6 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { None }; - // --- Input ------------------------------------------------------------- - let wanted = args.read_ids.as_deref().map(load_read_ids).transpose()?; - if let Some(w) = &wanted { - info!( - "aligning only the {} reads named in --read-ids", - style::count(w.len()) - ); - } - let threads = rayon::current_num_threads().max(1); - let (source, in_header) = match sniff(&args.reads)? { - InputFormat::Bam => { - let workers = std::num::NonZero::new(threads.div_ceil(4)).expect("at least one"); - let file = File::open(&args.reads)?; - let mut r = bam::io::Reader::from(bgzf::io::MultithreadedReader::with_worker_count( - workers, file, - )); - let h = r - .read_header() - .with_context(|| format!("reading BAM header of {}", args.reads.display()))?; - (Source::Bam(r), h) - } - InputFormat::Fasta => ( - Source::Fastx(FastxReader::new(open_text(&args.reads)?, false)), - sam::Header::default(), - ), - InputFormat::Fastq => ( - Source::Fastx(FastxReader::new(open_text(&args.reads)?, true)), - sam::Header::default(), - ), - }; let out_header = build_header(aligner.panel(), &in_header, args.sort)?; // --- Output ------------------------------------------------------------ @@ -1170,18 +1994,26 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { let batch_size = args.batch_size; let reader = std::thread::Builder::new() .name("eschalign-reader".into()) - .spawn(move || read_batches(source, wanted, batch_size, in_tx))?; + .spawn(move || match input { + Input::Single(source) => read_batches(source, wanted, batch_size, in_tx), + Input::Paired(source) => read_pair_batches(source, wanted, batch_size, in_tx), + })?; // With a GPU, its thread sits between the reader and the dispatcher; the // channel after it is as shallow as the one before, for the same memory // bound. let (batches, gpu_thread): (Box>, _) = match gpu_scorer { Some(scorer) => { let (tx, rx) = sync_channel::(2); - let both = opts.both_strands; + // A pair's fallback maps each mate on either strand, so a paired + // run always scores both orientations of every mate. + let both = opts.both_strands || paired; let max_read_len = args.max_read_len; + let gpu_pair = pair_cfg.clone(); let t = std::thread::Builder::new() .name("eschalign-gpu".into()) - .spawn(move || gpu::score_batches(scorer, in_rx, tx, both, max_read_len))?; + .spawn(move || { + gpu::score_batches(scorer, in_rx, tx, both, max_read_len, gpu_pair) + })?; (Box::new(rx.into_iter()), Some(t)) } None => (Box::new(in_rx.into_iter().map(|b| (b, None))), None), @@ -1204,6 +2036,9 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { split_strands: AtomicU64::new(0), max_read_len: args.max_read_len, skipped_long: AtomicU64::new(0), + pair: pair_cfg, + pairs: AtomicU64::new(0), + proper_pairs: AtomicU64::new(0), }; let busy = dispatch(&ctx, batches, res_tx, permit_rx); // Join every thread before deciding which error to report: a failed stage @@ -1238,7 +2073,12 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { ); } let summary = write_result.with_context(|| format!("writing {}", args.output.display()))?; - let counts = read_result.with_context(|| format!("reading {}", args.reads.display()))?; + let input_name = args + .reads + .as_deref() + .or(args.r1.as_deref()) + .map_or_else(String::new, |p| p.display().to_string()); + let counts = read_result.with_context(|| format!("reading {input_name}"))?; let skipped_long = ctx.skipped_long.load(Ordering::Relaxed); debug!( "workers were busy {:.1} s of {:.1} s x {} threads; {} read strands of {}+ nt \ @@ -1271,6 +2111,19 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { args.max_read_len ); } + if paired { + let (pairs, proper) = ( + ctx.pairs.load(Ordering::Relaxed), + ctx.proper_pairs.load(Ordering::Relaxed), + ); + info!( + "{} read pairs: {} proper ({:.1}%), {} not", + style::count(pairs), + style::count(proper), + 100.0 * proper as f64 / pairs.max(1) as f64, + style::count(pairs - proper), + ); + } let wall = started.elapsed().as_secs_f64(); let mapped = summary.unique + summary.tied; let pct = |n: u64| 100.0 * n as f64 / summary.reads.max(1) as f64; @@ -1295,21 +2148,36 @@ pub fn run(args: AlignArgs) -> anyhow::Result<()> { /// bases, whichever comes first — so a read of hundreds of kilobases (they /// occur: up to ~400 kb in a real tRNA sample) is a chunk of its own rather /// than the tail that sixty-three short reads wait behind. -fn cut_chunks(batch: Vec) -> Vec> { +/// +/// With `paired`, the batch is mate 1, mate 2, mate 1, … and the two mates of +/// a pair are one unit: a cut never falls between them, so a pair is always +/// aligned, and its records written, by one worker. +fn cut_chunks(batch: Vec, paired: bool) -> Vec> { + let len_of = |rec: &InRecord| match rec { + InRecord::Bam(r) => r.sequence().len(), + InRecord::Buf(r) => r.sequence().len(), + }; + let unit = if paired { 2 } else { 1 }; + debug_assert!( + batch.len().is_multiple_of(unit), + "a paired batch holds whole pairs" + ); let mut chunks = Vec::new(); let mut cur = Vec::with_capacity(CHUNK); let mut bases = 0usize; - for rec in batch { - let len = match &rec { - InRecord::Bam(r) => r.sequence().len(), - InRecord::Buf(r) => r.sequence().len(), - }; + let mut it = batch.into_iter(); + loop { + let group: Vec = it.by_ref().take(unit).collect(); + if group.is_empty() { + break; + } + let len: usize = group.iter().map(len_of).sum(); if !cur.is_empty() && (cur.len() >= CHUNK || bases + len > CHUNK_BASES) { chunks.push(std::mem::replace(&mut cur, Vec::with_capacity(CHUNK))); bases = 0; } bases += len; - cur.push(rec); + cur.extend(group); } if !cur.is_empty() { chunks.push(cur); @@ -1340,7 +2208,7 @@ fn dispatch( let mut seq = 0u64; 'batches: for (batch, scores) in batches { let mut first = 0usize; - for chunk in cut_chunks(batch) { + for chunk in cut_chunks(batch, ctx.pair.is_some()) { if permit_rx.recv().is_err() { // The writer stopped; its error is reported after the join. break 'batches; @@ -1352,7 +2220,10 @@ fn dispatch( let this = seq; scope.spawn(move |_| { let t0 = Instant::now(); - let r = process_chunk(ctx, chunk, chunk_scores); + let r = match &ctx.pair { + Some(cfg) => process_pair_chunk(ctx, cfg, chunk, chunk_scores), + None => process_chunk(ctx, chunk, chunk_scores), + }; busy_ns.fetch_add(t0.elapsed().as_nanos() as u64, Ordering::Relaxed); let _ = tx.send((this, r)); }); @@ -1823,6 +2694,7 @@ mod gpu { tx: SyncSender, both_strands: bool, max_read_len: usize, + pair: Option, ) -> anyhow::Result { use eschalign::alphabet::{encode_base, reverse_complement_codes}; let mut stats = GpuStats { @@ -1845,17 +2717,41 @@ mod gpu { // A read over --max-read-len goes as an empty query, which // the kernel leaves unscored; the worker never asks for it. let mut skipped = 0u64; - let mut codes: Vec> = batch - .iter() - .map(|rec| { - letters(rec, &mut buf); - if over_limit(max_read_len, buf.len()) { - skipped += 1; - return Vec::new(); + let mut codes: Vec> = match &pair { + None => batch + .iter() + .map(|rec| { + letters(rec, &mut buf); + if over_limit(max_read_len, buf.len()) { + skipped += 1; + return Vec::new(); + } + buf.iter().map(|&b| encode_base(b)).collect() + }) + .collect(), + // Pairs: the aligned bases are each mate less its 3' + // adapter, found by the same function the workers use. + Some(cfg) => { + let mut buf2 = Vec::new(); + let mut codes = Vec::with_capacity(batch.len()); + for mates in batch.chunks(2) { + letters(&mates[0], &mut buf); + letters(&mates[1], &mut buf2); + let ad = pair_adapters(&buf, &buf2, cfg, max_read_len); + for (l, ad) in [(&buf, ad[0]), (&buf2, ad[1])] { + if over_limit(max_read_len, l.len()) { + skipped += 1; + codes.push(Vec::new()); + } else { + codes.push( + l[..l.len() - ad].iter().map(|&b| encode_base(b)).collect(), + ); + } + } } - buf.iter().map(|&b| encode_base(b)).collect() - }) - .collect(); + codes + } + }; if both_strands { skipped *= 2; } @@ -1917,6 +2813,7 @@ mod gpu { _tx: SyncSender, _both_strands: bool, _max_read_len: usize, + _pair: Option, ) -> anyhow::Result { match scorer {} } @@ -1941,7 +2838,7 @@ mod tests { |len: usize| InRecord::Buf(fastx_to_record(b"r".to_vec(), vec![b'A'; len], None)); let mut batch: Vec = (0..100).map(|_| read(100)).collect(); batch.insert(10, read(300_000)); - let sizes: Vec = cut_chunks(batch).iter().map(Vec::len).collect(); + let sizes: Vec = cut_chunks(batch, false).iter().map(Vec::len).collect(); // 10 short, the giant alone, then the other 90 in chunks of <= CHUNK. assert_eq!(sizes, vec![10, 1, CHUNK, 90 - CHUNK]); } diff --git a/crates/eschalign-cli/tests/align_paired.rs b/crates/eschalign-cli/tests/align_paired.rs new file mode 100644 index 0000000..2037c60 --- /dev/null +++ b/crates/eschalign-cli/tests/align_paired.rs @@ -0,0 +1,805 @@ +// SPDX-License-Identifier: MIT + +//! Paired-end `eschalign` (rnabioco/eschalign#1) over the synthesized fixture +//! in `fixtures/pairs/` (regenerate with `gen_pair_fixture.py`): a four- +//! reference panel — two paralogs `P1`/`P2` that differ at one base, and two +//! unrelated references — and one read pair per behaviour, each named for what +//! it pins. No real data, no pod5. +//! +//! Every test aligns with `--min-score 40` unless it says otherwise: semi-global +//! scoring (the paired default) gives any read some positive score somewhere, so +//! a floor is what makes a mate that matches nothing "not aligned". + +use noodles_bam as bam; +use noodles_sam as sam; +use sam::alignment::RecordBuf; +use sam::alignment::io::Write as _; +use sam::alignment::record::Flags; +use sam::alignment::record::data::field::Tag; +use sam::alignment::record_buf::data::field::Value; +use std::path::{Path, PathBuf}; +use std::process::{Command, Output}; + +fn fx() -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures/pairs") +} + +fn panel() -> PathBuf { + fx().join("pairs_panel.fa") +} + +const STD: &[&str] = &["--min-score", "40"]; + +fn eschalign(args: &[&std::ffi::OsStr]) -> Output { + Command::new(env!("CARGO_BIN_EXE_eschalign")) + .args(args) + .output() + .expect("eschalign should launch") +} + +fn run_ok(args: &[&std::ffi::OsStr]) -> String { + let o = eschalign(args); + let stderr = String::from_utf8_lossy(&o.stderr).into_owned(); + assert!(o.status.success(), "eschalign failed:\n{stderr}"); + stderr +} + +/// `eschalign -1 R1 -2 R2 -r panel -o out `, returning the records. +fn align_files( + dir: &Path, + tag: &str, + r1: &Path, + r2: &Path, + extra: &[&str], +) -> (sam::Header, Vec) { + let out = dir.join(format!("{tag}.bam")); + let panel_path = panel(); + let mut args: Vec<&std::ffi::OsStr> = vec![ + "-1".as_ref(), + r1.as_os_str(), + "-2".as_ref(), + r2.as_os_str(), + "-r".as_ref(), + panel_path.as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]; + args.extend(extra.iter().map(|s| std::ffi::OsStr::new(*s))); + run_ok(&args); + read_bam(&out) +} + +fn align_fixture(dir: &Path, tag: &str, extra: &[&str]) -> (sam::Header, Vec) { + let mut all = STD.to_vec(); + all.extend_from_slice(extra); + align_files( + dir, + tag, + &fx().join("pairs_R1.fq"), + &fx().join("pairs_R2.fq"), + &all, + ) +} + +fn read_bam(path: &Path) -> (sam::Header, Vec) { + let mut r = bam::io::Reader::new(std::fs::File::open(path).unwrap()); + let header = r.read_header().unwrap(); + let recs = r + .record_bufs(&header) + .collect::, _>>() + .unwrap(); + (header, recs) +} + +fn name(r: &RecordBuf) -> String { + r.name().map(|n| n.to_string()).unwrap_or_default() +} + +fn ref_name(h: &sam::Header, r: &RecordBuf) -> Option { + r.reference_sequence_id() + .map(|i| h.reference_sequences().get_index(i).unwrap().0.to_string()) +} + +fn cigar(r: &RecordBuf) -> String { + use sam::alignment::record::cigar::op::Kind; + r.cigar() + .as_ref() + .iter() + .map(|op| { + let c = match op.kind() { + Kind::Match => 'M', + Kind::Insertion => 'I', + Kind::Deletion => 'D', + Kind::SoftClip => 'S', + _ => '?', + }; + format!("{}{c}", op.len()) + }) + .collect() +} + +fn string_tag(r: &RecordBuf, tag: &[u8; 2]) -> Option { + match r.data().get(&Tag::from(*tag)) { + Some(Value::String(s)) => Some(s.to_string()), + _ => None, + } +} + +fn int_tag(r: &RecordBuf, tag: &[u8; 2]) -> Option { + r.data().get(&Tag::from(*tag)).and_then(|v| v.as_int()) +} + +fn flag(r: &RecordBuf) -> u16 { + u16::from(r.flags()) +} + +fn pos(r: &RecordBuf) -> Option { + r.alignment_start().map(|p| p.get()) +} + +/// The two records of the pair called `pair`, mate 1 then mate 2, which must +/// be adjacent and in that order. +fn mates<'a>(recs: &'a [RecordBuf], pair: &str) -> [&'a RecordBuf; 2] { + let i = recs + .iter() + .position(|r| name(r) == pair && r.flags().is_first_segment()) + .unwrap_or_else(|| panic!("no mate 1 named {pair}")); + let (a, b) = (&recs[i], &recs[i + 1]); + assert_eq!(name(b), pair, "mates are adjacent"); + assert!(b.flags().is_last_segment(), "mate 2 follows mate 1"); + [a, b] +} + +fn fasta_seq(refname: &str) -> String { + let text = std::fs::read_to_string(panel()).unwrap(); + let mut it = text.lines(); + while let Some(l) = it.next() { + if l == format!(">{refname}") { + return it.next().unwrap().to_string(); + } + } + panic!("no reference {refname}"); +} + +/// The two mates' alignments, flags and mate fields, as one comparable line +/// each (everything but the `@PG`). +fn lines(recs: &[RecordBuf]) -> Vec { + recs.iter() + .map(|r| { + let mut tags: Vec = r + .data() + .iter() + .map(|(t, v)| format!("{}{}={v:?}", t.as_ref()[0] as char, t.as_ref()[1] as char)) + .collect(); + tags.sort(); + format!( + "{} {} {:?} {:?} {:?} {} {:?} {:?} {} {} {} {}", + name(r), + flag(r), + r.reference_sequence_id(), + r.alignment_start(), + r.mapping_quality(), + cigar(r), + r.mate_reference_sequence_id(), + r.mate_alignment_start(), + r.template_length(), + String::from_utf8_lossy(r.sequence().as_ref()), + r.quality_scores() + .as_ref() + .iter() + .map(|q| (q + 33) as char) + .collect::(), + tags.join(",") + ) + }) + .collect() +} + +#[test] +fn pair_resolves_paralog_tie() { + let dir = tempfile::tempdir().unwrap(); + let (h, recs) = align_fixture(dir.path(), "paired", &[]); + let [a, b] = mates(&recs, "paralog"); + for r in [a, b] { + assert_eq!(ref_name(&h, r).as_deref(), Some("P1"), "{}", name(r)); + assert_eq!(r.mapping_quality().map(|q| q.get()), Some(60)); + assert!(string_tag(r, b"XA").is_none(), "a unique pair has no XA"); + } + // Mate 1 alone is a tie between the paralogs: the unpaired aligner says so. + let out = dir.path().join("single.bam"); + run_ok(&[ + fx().join("pairs_R1.fq").as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + "--min-score".as_ref(), + "40".as_ref(), + ]); + let (_, single) = read_bam(&out); + let alone = single.iter().find(|r| name(r) == "paralog").unwrap(); + assert_eq!(alone.mapping_quality().map(|q| q.get()), Some(0)); + assert!(string_tag(alone, b"XA").is_some()); + // The pair's mate 1 still reports the tie it would have had alone as XS. + assert_eq!(int_tag(a, b"XS"), int_tag(a, b"AS")); + // Flags: paired, proper, mate 1 forward / mate 2 reverse. + assert_eq!(flag(a), 0x1 | 0x2 | 0x20 | 0x40); + assert_eq!(flag(b), 0x1 | 0x2 | 0x10 | 0x80); + assert_eq!((pos(a), pos(b)), (Some(21), Some(91))); + assert_eq!( + ( + a.mate_alignment_start().map(|p| p.get()), + a.template_length() + ), + (Some(91), 130) + ); + assert_eq!(b.template_length(), -130); + assert_eq!(string_tag(a, b"MC").as_deref(), Some("60M")); + // Mate 2 is stored reverse-complemented with its qualities reversed. + let (_, input) = { + let text = std::fs::read_to_string(fx().join("pairs_R2.fq")).unwrap(); + let l: Vec<&str> = text.lines().collect(); + ((), (l[1].to_string(), l[3].to_string())) + }; + let rev_q: String = input.1.chars().rev().collect(); + let got_q: String = b + .quality_scores() + .as_ref() + .iter() + .map(|q| (q + 33) as char) + .collect(); + assert_eq!(got_q, rev_q); +} + +#[test] +fn pair_tie_is_reported_in_xa_and_secondaries() { + let dir = tempfile::tempdir().unwrap(); + let (h, recs) = align_fixture(dir.path(), "tie", &[]); + let [a, b] = mates(&recs, "pair_tie"); + for r in [a, b] { + assert_eq!(r.mapping_quality().map(|q| q.get()), Some(0)); + assert_eq!(ref_name(&h, r).as_deref(), Some("P1")); + let xa = string_tag(r, b"XA").expect("XA"); + assert!(xa.starts_with("P2,"), "{xa}"); + } + // The k-th XA entry of each mate belongs to the same other pair. + let (_, sec) = align_fixture(dir.path(), "tie_sec", &["--secondary"]); + let secs: Vec<&RecordBuf> = sec + .iter() + .filter(|r| name(r) == "pair_tie" && r.flags().is_secondary()) + .collect(); + assert_eq!(secs.len(), 2, "one secondary per mate per further pair"); + assert!(secs[0].flags().is_first_segment() && secs[1].flags().is_last_segment()); + assert_eq!(ref_name(&h, secs[0]).as_deref(), Some("P2")); + assert_eq!( + secs[0].mate_reference_sequence_id(), + secs[1].reference_sequence_id() + ); + assert!(secs[0].flags().is_properly_segmented()); +} + +#[test] +fn pair_insert_shorter_than_read() { + let dir = tempfile::tempdir().unwrap(); + let (h, recs) = align_fixture(dir.path(), "short", &[]); + let [a, b] = mates(&recs, "short_insert"); + assert_eq!(ref_name(&h, a).as_deref(), Some("U1")); + // Insert 40, reads 60: TLEN is the insert, and each read's 20 adapter + // bases are a soft clip on its 3' end (the start of the reverse mate). + assert_eq!((a.template_length(), b.template_length()), (40, -40)); + assert_eq!(cigar(a), "40M20S"); + assert_eq!(cigar(b), "20S40M"); + assert_eq!((pos(a), pos(b)), (Some(31), Some(31))); + assert!(a.flags().is_properly_segmented()); + // MC is the *mate's* CIGAR. + assert_eq!(string_tag(a, b"MC").as_deref(), Some("20S40M")); + assert_eq!(string_tag(b, b"MC").as_deref(), Some("40M20S")); + // Nothing discarded: SEQ keeps the adapter. + assert_eq!(a.sequence().len(), 60); + assert_eq!(b.sequence().len(), 60); +} + +#[test] +fn pair_adapter_tag_counts() { + let dir = tempfile::tempdir().unwrap(); + let (_, recs) = align_fixture(dir.path(), "ad", &[]); + let ad = |pair: &str| { + let [a, b] = mates(&recs, pair); + (int_tag(a, b"ad"), int_tag(b, b"ad")) + }; + assert_eq!(ad("short_insert"), (Some(20), Some(20))); + // One mate reads through: its own count, and 0 (present) on the other. + assert_eq!(ad("one_mate_reads_through"), (Some(15), Some(0))); + assert_eq!(ad("plain"), (Some(0), Some(0))); + let [a, b] = mates(&recs, "one_mate_reads_through"); + assert_eq!(cigar(a), "40M15S"); + assert_eq!(cigar(b), "35M"); + assert_eq!(a.template_length(), 40); + // --no-adapter-trim: no detection, no tag, and the adapter is aligned (or + // overhangs) rather than clipped by us. + let (_, off) = align_fixture(dir.path(), "noad", &["--no-adapter-trim"]); + assert!(off.iter().all(|r| int_tag(r, b"ad").is_none())); + // Custom adapters: an adapter that is not in the data finds nothing. + let (_, custom) = align_fixture( + dir.path(), + "custom", + &[ + "--adapter-r1", + "TTTTTTTTTTTTTTTT", + "--adapter-r2", + "TTTTTTTTTTTTTTTT", + ], + ); + let [a, _] = mates(&custom, "short_insert"); + assert_eq!(int_tag(a, b"ad"), Some(0)); + // A raised minimum overlap hides the 15-base read-through. + let (_, hi) = align_fixture(dir.path(), "minov", &["--min-adapter-overlap", "16"]); + let [a, _] = mates(&hi, "one_mate_reads_through"); + assert_eq!(int_tag(a, b"ad"), Some(0)); +} + +#[test] +fn pair_terminal_mismatch_not_clipped() { + let dir = tempfile::tempdir().unwrap(); + let (_, recs) = align_fixture(dir.path(), "tm", &[]); + let [a, b] = mates(&recs, "terminal_mismatch"); + let u1 = fasta_seq("U1"); + // Mate 1: last base is a mismatch, kept — MD ends in it, NM counts it. + assert_eq!(cigar(a), "60M"); + assert_eq!(int_tag(a, b"NM"), Some(1)); + let last = u1.as_bytes()[69] as char; + assert_eq!(string_tag(a, b"MD"), Some(format!("59{last}0"))); + // Mate 2: its first base (the *last* of the reverse-complemented + // alignment) is a mismatch. + assert_eq!(cigar(b), "60M"); + assert_eq!(int_tag(b, b"NM"), Some(1)); + assert!(string_tag(b, b"MD").unwrap().ends_with('0')); + // The unpaired default (--mode local) would clip it; the paired default + // keeps it, and an explicit --mode local gives the clip back. + let (_, local) = align_fixture(dir.path(), "tm_local", &["--mode", "local"]); + let [a, _] = mates(&local, "terminal_mismatch"); + assert_eq!(cigar(a), "59M1S"); +} + +#[test] +fn pair_over_max_insert_not_proper() { + let dir = tempfile::tempdir().unwrap(); + let (h, recs) = align_fixture(dir.path(), "wide", &["--max-insert", "120"]); + let [a, b] = mates(&recs, "over_max_insert"); + assert!(!a.flags().is_properly_segmented() && !b.flags().is_properly_segmented()); + // Each mate at its own best reference, with the mate fields still set. + assert_eq!(ref_name(&h, a).as_deref(), Some("U1")); + assert_eq!(ref_name(&h, b).as_deref(), Some("U1")); + assert_eq!(flag(a), 0x1 | 0x20 | 0x40); + assert_eq!(flag(b), 0x1 | 0x10 | 0x80); + assert_eq!(a.template_length(), 150); + // The default --max-insert is 500: the same pair is proper. + let (_, recs) = align_fixture(dir.path(), "roomy", &[]); + let [a, b] = mates(&recs, "over_max_insert"); + assert!(a.flags().is_properly_segmented() && b.flags().is_properly_segmented()); +} + +#[test] +fn pair_wrong_orientation_not_proper() { + let dir = tempfile::tempdir().unwrap(); + let (h, recs) = align_fixture(dir.path(), "ff", &[]); + let [a, b] = mates(&recs, "wrong_orientation"); + assert!(!a.flags().is_properly_segmented() && !b.flags().is_properly_segmented()); + // Both mates sit where they are, forward, at their own best reference. + assert_eq!(ref_name(&h, a).as_deref(), Some("U2")); + assert_eq!(ref_name(&h, b).as_deref(), Some("U2")); + assert_eq!(flag(a), 0x1 | 0x40); + assert_eq!(flag(b), 0x1 | 0x80); + // A mate that matches nothing is unmapped and rides at its partner's place. + let [a, b] = mates(&recs, "one_unmapped"); + assert_eq!(flag(a), 0x1 | 0x8 | 0x40); + assert_eq!(flag(b), 0x1 | 0x4 | 0x80); + assert_eq!(b.reference_sequence_id(), a.reference_sequence_id()); + assert_eq!(pos(b), pos(a)); + assert_eq!(b.mapping_quality().map(|q| q.get()), Some(0)); +} + +fn fastq_records(path: &Path) -> Vec<(String, String, String)> { + let text = std::fs::read_to_string(path).unwrap(); + let l: Vec<&str> = text.lines().collect(); + l.chunks(4) + .map(|c| (c[0][1..].to_string(), c[1].to_string(), c[3].to_string())) + .collect() +} + +/// The fixture as an interleaved FASTQ and an unaligned BAM. +fn other_forms(dir: &Path) -> (PathBuf, PathBuf) { + let (r1, r2) = ( + fastq_records(&fx().join("pairs_R1.fq")), + fastq_records(&fx().join("pairs_R2.fq")), + ); + let il = dir.join("pairs_interleaved.fq"); + let mut text = String::new(); + let ubam = dir.join("pairs.ubam"); + let header = sam::Header::default(); + let mut w = bam::io::Writer::new(std::fs::File::create(&ubam).unwrap()); + w.write_header(&header).unwrap(); + for (a, b) in r1.iter().zip(&r2) { + for (m, (n, s, q)) in [a, b].into_iter().enumerate() { + text.push_str(&format!("@{n}\n{s}\n+\n{q}\n")); + let mut flags = Flags::SEGMENTED | Flags::UNMAPPED | Flags::MATE_UNMAPPED; + flags |= if m == 0 { + Flags::FIRST_SEGMENT + } else { + Flags::LAST_SEGMENT + }; + let rec = RecordBuf::builder() + .set_name(n.clone().into_bytes()) + .set_flags(flags) + .set_sequence(s.clone().into_bytes().into()) + .set_quality_scores(q.bytes().map(|c| c - 33).collect::>().into()) + .build(); + w.write_alignment_record(&header, &rec).unwrap(); + } + } + w.try_finish().unwrap(); + std::fs::write(&il, text).unwrap(); + (il, ubam) +} + +fn align_single_input(dir: &Path, tag: &str, reads: &Path, extra: &[&str]) -> Vec { + let out = dir.join(format!("{tag}.bam")); + let panel_path = panel(); + let mut args: Vec<&std::ffi::OsStr> = vec![ + reads.as_os_str(), + "-r".as_ref(), + panel_path.as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]; + args.extend(STD.iter().map(|s| std::ffi::OsStr::new(*s))); + args.extend(extra.iter().map(|s| std::ffi::OsStr::new(*s))); + run_ok(&args); + read_bam(&out).1 +} + +#[test] +fn pair_inputs_agree_fastq_interleaved_ubam() { + let dir = tempfile::tempdir().unwrap(); + let (il, ubam) = other_forms(dir.path()); + for extra in [&[][..], &["--strand", "both", "--secondary"][..]] { + let (_, files) = align_fixture(dir.path(), "files", extra); + let inter = + align_single_input(dir.path(), "il", &il, &[&["--interleaved"], extra].concat()); + let from_bam = align_single_input(dir.path(), "ubam", &ubam, extra); + assert!(files.len() >= 18); + assert_eq!( + lines(&files), + lines(&inter), + "interleaved differs {extra:?}" + ); + assert_eq!(lines(&files), lines(&from_bam), "uBAM differs {extra:?}"); + } + // Gzipped mate files, with /1 and /2 on the names, give the same records. + let gz = |src: &Path, suffix: &str, to: &Path| { + let mut text = String::new(); + for (n, s, q) in fastq_records(src) { + text.push_str(&format!("@{n}/{suffix} extra words\n{s}\n+\n{q}\n")); + } + let mut e = flate2::write::GzEncoder::new( + std::fs::File::create(to).unwrap(), + flate2::Compression::fast(), + ); + std::io::Write::write_all(&mut e, text.as_bytes()).unwrap(); + e.finish().unwrap(); + }; + let (g1, g2) = (dir.path().join("a_1.fq.gz"), dir.path().join("a_2.fq.gz")); + gz(&fx().join("pairs_R1.fq"), "1", &g1); + gz(&fx().join("pairs_R2.fq"), "2", &g2); + let mut all = STD.to_vec(); + all.push("--no-adapter-trim"); + let (_, suffixed) = align_files(dir.path(), "suffixed", &g1, &g2, &all); + let (_, plain) = align_fixture(dir.path(), "plain", &["--no-adapter-trim"]); + assert_eq!(lines(&suffixed), lines(&plain)); +} + +#[test] +fn pair_orphan_is_error() { + let dir = tempfile::tempdir().unwrap(); + let r1 = fastq_records(&fx().join("pairs_R1.fq")); + let r2 = fastq_records(&fx().join("pairs_R2.fq")); + let write = |p: &Path, recs: &[(String, String, String)]| { + let t: String = recs + .iter() + .map(|(n, s, q)| format!("@{n}\n{s}\n+\n{q}\n")) + .collect(); + std::fs::write(p, t).unwrap(); + }; + let fail = |args: &[&std::ffi::OsStr]| -> String { + let o = eschalign(args); + assert!(!o.status.success(), "an orphan mate must fail the run"); + String::from_utf8_lossy(&o.stderr).into_owned() + }; + let out = dir.path().join("o.bam"); + // Mate file one record short: the read left over is named. + let (short, full) = (dir.path().join("short.fq"), dir.path().join("full.fq")); + write(&short, &r2[..r2.len() - 1]); + write(&full, &r1); + let err = fail(&[ + "-1".as_ref(), + full.as_os_str(), + "-2".as_ref(), + short.as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]); + assert!(err.contains("plain") && err.contains("no mate"), "{err}"); + // Names that do not match at the same record. + let mut swapped = r2.clone(); + swapped.swap(0, 1); + let sw = dir.path().join("swapped.fq"); + write(&sw, &swapped); + let err = fail(&[ + "-1".as_ref(), + full.as_os_str(), + "-2".as_ref(), + sw.as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]); + assert!(err.contains("paralog") && err.contains("pair_tie"), "{err}"); + // Interleaved with a trailing orphan. + let il = dir.path().join("il.fq"); + let mut recs: Vec<_> = r1 + .iter() + .zip(&r2) + .flat_map(|(a, b)| [a.clone(), b.clone()]) + .collect(); + recs.push(( + "lonely".into(), + "ACGTACGTACGTACGTACGT".into(), + "I".repeat(20), + )); + write(&il, &recs); + let err = fail(&[ + il.as_os_str(), + "--interleaved".as_ref(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]); + assert!(err.contains("lonely") && err.contains("no mate"), "{err}"); + // A paired BAM whose last record has lost its mate. + let (_, ubam) = other_forms(dir.path()); + let (header, mut recs) = read_bam(&ubam); + recs.pop(); + let lost = dir.path().join("lost.bam"); + let mut w = bam::io::Writer::new(std::fs::File::create(&lost).unwrap()); + w.write_header(&header).unwrap(); + for r in &recs { + w.write_alignment_record(&header, r).unwrap(); + } + w.try_finish().unwrap(); + let err = fail(&[ + lost.as_os_str(), + "-r".as_ref(), + panel().as_os_str(), + "-o".as_ref(), + out.as_os_str(), + ]); + assert!(err.contains("plain") && err.contains("no mate"), "{err}"); +} + +#[test] +fn pair_ubam_tags_are_copied_on_both_mates() { + let dir = tempfile::tempdir().unwrap(); + let (_, ubam) = other_forms(dir.path()); + let (header, mut recs) = read_bam(&ubam); + for (k, r) in recs.iter_mut().enumerate() { + r.data_mut().insert( + Tag::new(b'Z', b'Z'), + Value::String(format!("tag{k}").into()), + ); + // A stale tag the aligner owns must be replaced, not kept. + r.data_mut().insert(Tag::new(b'N', b'M'), Value::UInt8(99)); + } + let tagged = dir.path().join("tagged.bam"); + let mut w = bam::io::Writer::new(std::fs::File::create(&tagged).unwrap()); + w.write_header(&header).unwrap(); + for r in &recs { + w.write_alignment_record(&header, r).unwrap(); + } + w.try_finish().unwrap(); + let out = align_single_input(dir.path(), "tagged_out", &tagged, &[]); + assert_eq!(out.len(), recs.len()); + for (k, r) in out.iter().enumerate() { + assert_eq!(string_tag(r, b"ZZ"), Some(format!("tag{k}"))); + assert_ne!(int_tag(r, b"NM"), Some(99)); + } +} + +#[test] +fn pair_sort_coordinate_matches_samtools() { + let dir = tempfile::tempdir().unwrap(); + for extra in [&[][..], &["--secondary", "--strand", "both"][..]] { + let (_, unsorted) = align_fixture(dir.path(), "u", extra); + let mut sorted_args = extra.to_vec(); + sorted_args.extend(["--sort", "coordinate"]); + let (sh, sorted) = align_fixture(dir.path(), "s", &sorted_args); + // The model: a stable sort on (tid, pos, reverse), no reference last. + let mut model = unsorted.clone(); + model.sort_by_key(|r| { + ( + r.reference_sequence_id().unwrap_or(usize::MAX), + r.alignment_start().map_or(0, |p| p.get()), + r.flags().is_reverse_complemented(), + ) + }); + assert_eq!(lines(&sorted), lines(&model), "model order {extra:?}"); + // Mates of a pair stay consistent: each record's mate fields point at + // its partner's place, wherever the sort put them. + for r in sorted.iter().filter(|r| !r.flags().is_secondary()) { + let partner = sorted + .iter() + .find(|p| { + name(p) == name(r) + && !p.flags().is_secondary() + && p.flags().is_first_segment() != r.flags().is_first_segment() + }) + .expect("partner"); + if !r.flags().is_mate_unmapped() { + assert_eq!(r.mate_alignment_start(), partner.alignment_start()); + assert_eq!( + r.mate_reference_sequence_id(), + partner.reference_sequence_id() + ); + } + } + assert!(sh.header().is_some()); + // And samtools itself, when there is one. + let Some(st) = samtools() else { + eprintln!( + "[align_paired] pair_sort_coordinate_matches_samtools: no samtools; model only" + ); + continue; + }; + let theirs = dir.path().join("theirs.bam"); + let o = Command::new(st) + .args(["sort", "--no-PG", "-o"]) + .arg(&theirs) + .arg(dir.path().join("u.bam")) + .output() + .unwrap(); + assert!(o.status.success(), "{}", String::from_utf8_lossy(&o.stderr)); + let (_, st_recs) = read_bam(&theirs); + assert_eq!(lines(&sorted), lines(&st_recs), "samtools order {extra:?}"); + } +} + +// SAMTOOLS and PATH are third-party variables, not ESCHALIGN_* knobs. +#[allow(clippy::disallowed_methods)] +fn samtools() -> Option { + if let Some(p) = std::env::var_os("SAMTOOLS") { + return Some(PathBuf::from(p)); + } + std::env::split_paths(&std::env::var_os("PATH")?) + .map(|d| d.join("samtools")) + .find(|p| p.is_file()) +} + +/// The GPU and CPU scorers feed the same pair resolution, tracebacks and +/// record building, so the records must be identical. +#[cfg(feature = "gpu")] +#[test] +fn pair_gpu_matches_cpu() { + let dir = tempfile::tempdir().unwrap(); + for (k, extra) in [ + &[][..], + &["--strand", "both", "--secondary"][..], + &["--max-insert", "120", "--batch-size", "4"][..], + &["--no-adapter-trim", "--mode", "local"][..], + ] + .iter() + .enumerate() + { + let (r1, r2, panel_path) = (fx().join("pairs_R1.fq"), fx().join("pairs_R2.fq"), panel()); + let mk = |device: &str, tag: &str| { + let out = dir.path().join(format!("{tag}{k}.bam")); + let mut args: Vec<&std::ffi::OsStr> = vec![ + "-1".as_ref(), + r1.as_os_str(), + "-2".as_ref(), + r2.as_os_str(), + "-r".as_ref(), + panel_path.as_os_str(), + "-o".as_ref(), + out.as_os_str(), + "--device".as_ref(), + device.as_ref(), + ]; + args.extend(STD.iter().map(|s| std::ffi::OsStr::new(*s))); + args.extend(extra.iter().map(|s| std::ffi::OsStr::new(*s))); + (eschalign(&args), out.clone()) + }; + let (g, gout) = mk("gpu", "gpu"); + let stderr = String::from_utf8_lossy(&g.stderr); + if !g.status.success() { + assert!( + stderr.contains("--device gpu cannot run"), + "GPU run failed for a reason other than a missing device:\n{stderr}" + ); + eprintln!("[align_paired] skipping pair_gpu_matches_cpu: {stderr}"); + return; + } + assert!(stderr.contains("on GPU"), "{stderr}"); + let (c, cout) = mk("cpu", "cpu"); + assert!(c.status.success(), "{}", String::from_utf8_lossy(&c.stderr)); + let (_, cpu) = read_bam(&cout); + let (_, gpu) = read_bam(&gout); + assert!(!cpu.is_empty()); + assert_eq!(cpu, gpu, "flags {extra:?}"); + } +} + +/// Chunks and batches are cut between pairs, never inside one: a few hundred +/// pairs through tiny batches on several threads come out in input order, mates +/// adjacent, record for record what one big batch gives. +#[test] +fn pair_chunking_keeps_mates_together() { + let dir = tempfile::tempdir().unwrap(); + let refs: Vec = ["P1", "P2", "U1", "U2"].map(fasta_seq).to_vec(); + let rc = |s: &str| -> String { + s.bytes() + .rev() + .map(|b| match b { + b'A' => 'T', + b'C' => 'G', + b'G' => 'C', + _ => 'A', + }) + .collect() + }; + let mut x = 0x9e37_79b9_7f4a_7c15u64; + let mut next = |n: usize| { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + (x % n as u64) as usize + }; + let (mut t1, mut t2) = (String::new(), String::new()); + for k in 0..300 { + let r = &refs[next(refs.len())]; + let start = next(40); + let len = 30 + next(40); + let insert = &r[start..start + len + next(30)]; + let q = "I".repeat(len); + t1.push_str(&format!("@p{k}\n{}\n+\n{q}\n", &insert[..len])); + t2.push_str(&format!( + "@p{k}\n{}\n+\n{q}\n", + rc(&insert[insert.len() - len..]) + )); + } + let (f1, f2) = (dir.path().join("r1.fq"), dir.path().join("r2.fq")); + std::fs::write(&f1, t1).unwrap(); + std::fs::write(&f2, t2).unwrap(); + let (_, big) = align_files(dir.path(), "big", &f1, &f2, STD); + let (_, small) = align_files( + dir.path(), + "small", + &f1, + &f2, + &[STD, &["--batch-size", "7", "-t", "4", "--secondary"]].concat(), + ); + assert_eq!(big.len(), 600); + let primary: Vec = small + .iter() + .filter(|r| !r.flags().is_secondary()) + .cloned() + .collect(); + assert_eq!(lines(&big), lines(&primary)); + for (k, pair) in big.chunks(2).enumerate() { + assert_eq!(name(&pair[0]), format!("p{k}")); + assert_eq!(name(&pair[1]), format!("p{k}")); + assert!(pair[0].flags().is_first_segment() && pair[1].flags().is_last_segment()); + } +} diff --git a/crates/eschalign-cli/tests/align_unpaired.rs b/crates/eschalign-cli/tests/align_unpaired.rs new file mode 100644 index 0000000..c171bf0 --- /dev/null +++ b/crates/eschalign-cli/tests/align_unpaired.rs @@ -0,0 +1,100 @@ +// SPDX-License-Identifier: MIT + +//! `unpaired_output_unchanged`: the unpaired path must write the same bytes it +//! wrote before paired-end support existed. +//! +//! The golden (`fixtures/unpaired_golden.tsv`) was measured with the binary at +//! the commit before `feat/paired-end`: for each flag combination, the record +//! count and an FNV-1a 64 digest of the *uncompressed* BAM record stream (the +//! header is skipped — its `@PG CL` carries the temporary paths of the run). +//! `ESCHALIGN_BLESS=1` rewrites it; do that only when an unpaired change is +//! intended, and say so in the commit. + +use std::io::Read; +use std::path::{Path, PathBuf}; +use std::process::Command; + +fn fixtures() -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/fixtures") +} + +const CONFIGS: &[(&str, &[&str])] = &[ + ("default", &[]), + ( + "semiglobal_both_secondary", + &["--mode", "semiglobal", "--strand", "both", "--secondary"], + ), + ("sort_coordinate", &["--sort", "coordinate"]), + ( + "min_score_max_ties", + &["--min-score", "60", "--max-ties", "1"], + ), + ( + "scoring_bwa", + &["--scoring", "1,-1,-2,-1", "--strand", "both"], + ), + ("max_read_len", &["--max-read-len", "60"]), +]; + +/// `(records, digest)` of a BAM's record stream, header skipped. +fn digest(path: &Path) -> (usize, u64) { + let mut raw = Vec::new(); + noodles_bgzf::io::Reader::new(std::fs::File::open(path).unwrap()) + .read_to_end(&mut raw) + .unwrap(); + let u32_at = |o: usize| u32::from_le_bytes(raw[o..o + 4].try_into().unwrap()) as usize; + assert_eq!(&raw[..4], b"BAM\x01"); + let mut off = 8 + u32_at(4); + let n_ref = u32_at(off); + off += 4; + for _ in 0..n_ref { + off += 4 + u32_at(off) + 4; + } + let body = &raw[off..]; + let (mut n, mut p) = (0, 0); + while p < body.len() { + p += 4 + u32::from_le_bytes(body[p..p + 4].try_into().unwrap()) as usize; + n += 1; + } + let mut h = 0xcbf2_9ce4_8422_2325u64; + for &b in body { + h ^= u64::from(b); + h = h.wrapping_mul(0x0000_0100_0000_01b3); + } + (n, h) +} + +#[test] +fn unpaired_output_unchanged() { + let dir = tempfile::tempdir().unwrap(); + let mut got = String::new(); + for (label, extra) in CONFIGS { + let out = dir.path().join(format!("{label}.bam")); + let o = Command::new(env!("CARGO_BIN_EXE_eschalign")) + .arg(fixtures().join("trna_mappings_padded.bam")) + .arg("-r") + .arg(fixtures().join("trna_reference.fa")) + .arg("-o") + .arg(&out) + .args(*extra) + .output() + .expect("eschalign should launch"); + assert!( + o.status.success(), + "{label}: {}", + String::from_utf8_lossy(&o.stderr) + ); + let (n, h) = digest(&out); + got.push_str(&format!("{label}\t{n}\t{h:016x}\n")); + } + let golden = fixtures().join("unpaired_golden.tsv"); + // A test-only switch, not a runtime knob: any value blesses. + #[allow(clippy::disallowed_methods)] + let bless = std::env::var_os("ESCHALIGN_BLESS").is_some(); + if bless { + std::fs::write(&golden, &got).unwrap(); + return; + } + let want = std::fs::read_to_string(&golden).expect("fixtures/unpaired_golden.tsv"); + assert_eq!(got, want, "the unpaired output changed"); +} diff --git a/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py b/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py new file mode 100644 index 0000000..debe94d --- /dev/null +++ b/crates/eschalign-cli/tests/fixtures/pairs/gen_pair_fixture.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +"""Regenerate the paired-end fixture: pairs_panel.fa, pairs_R1.fq, pairs_R2.fq. + +Everything is synthesized from a fixed seed (no real data). The truth set is +built by hand, one pair per behaviour the paired-end tests pin; `align_paired.rs` +states what each pair must produce. Run from this directory: + + python3 gen_pair_fixture.py +""" +import random + +R1_ADAPTER = "AGATCGGAAGAGCACACGTCTGAACTCCAGTCA" +R2_ADAPTER = "AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT" + +rng = random.Random(20261004) + + +def rand(n): + return "".join(rng.choice("ACGT") for _ in range(n)) + + +def rc(s): + return s[::-1].translate(str.maketrans("ACGT", "TGCA")) + + +def other(b): + return {"A": "C", "C": "G", "G": "T", "T": "A"}[b] + + +# Two paralogs differing at one base (position 100), and two unrelated references. +p1 = rand(150) +p2 = p1[:100] + other(p1[100]) + p1[101:] +u1 = rand(150) +u2 = rand(150) +panel = {"P1": p1, "P2": p2, "U1": u1, "U2": u2} + +pairs = [] + + +def add(name, r1, r2): + pairs.append((name, r1, r2)) + + +# Mate 1 lies where P1 and P2 are identical (alone it ties); mate 2 covers +# position 100, which only P1 explains. +add("paralog", p1[20:80], rc(p1[90:150])) +# Both mates in the identical region: the pair ties between P1 and P2. +add("pair_tie", p1[0:50], rc(p1[40:90])) +# Insert of 40 nt, both mates read 20 nt into their adapter. +tpl = u1[30:70] +add("short_insert", tpl + R1_ADAPTER[:20], rc(tpl) + R2_ADAPTER[:20]) +# Only mate 1 reads through (15 adapter bases); mate 2 is shorter and stays inside. +add("one_mate_reads_through", u2[30:70] + R1_ADAPTER[:15], rc(u2[35:70])) +# A mismatch on the last base of mate 1 and on the first base of mate 2. +m1 = u1[10:70] +m1 = m1[:-1] + other(m1[-1]) +m2 = rc(u1[70:130]) +m2 = other(m2[0]) + m2[1:] +add("terminal_mismatch", m1, m2) +# Insert of 150: proper by default, not under --max-insert 120. +add("over_max_insert", u1[0:50], rc(u1[100:150])) +# Both mates read as sequenced on the same strand (FF). +add("wrong_orientation", u2[0:50], u2[100:150]) +# Mate 2 is junk. +add("one_unmapped", u2[20:80], rand(60)) +# An ordinary pair. +add("plain", u2[10:70], rc(u2[80:140])) + + +def qual(n, shift): + return "".join(chr(33 + 20 + (i * 3 + shift) % 20) for i in range(n)) + + +with open("pairs_panel.fa", "w") as f: + for k, v in panel.items(): + f.write(f">{k}\n{v}\n") +with open("pairs_R1.fq", "w") as f1, open("pairs_R2.fq", "w") as f2: + for name, a, b in pairs: + f1.write(f"@{name}\n{a}\n+\n{qual(len(a), 0)}\n") + f2.write(f"@{name}\n{b}\n+\n{qual(len(b), 7)}\n") diff --git a/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq b/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq new file mode 100644 index 0000000..8cd2326 --- /dev/null +++ b/crates/eschalign-cli/tests/fixtures/pairs/pairs_R1.fq @@ -0,0 +1,36 @@ +@paralog +CCACTTCGACGTGGGGAAGCCGTACCGGTGACAGGGTTTGTCGTAATGTGACCGGGGTGC ++ +58;>ADG69ADG69ADG69ADG69ADG69ADG69< +@short_insert +ATGAGTGAGCCATCCAACATACGTGGCCGGCCGGGTCTAGAGATCGGAAGAGCACACGTC ++ +58;>ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69< +@wrong_orientation +TTTTGTCGACAAACGTCCGTCTTGCGTATGGCACCTTTCCGGTATTGCTA ++ +58;>ADG69ADG69ADG69< +@one_unmapped +CTTGCGTATGGCACCTTTCCGGTATTGCTACAACCATCTGGGGGGCCGGATCTATTCTGG ++ +58;>ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69ADG69 +@pair_tie +GATGGAAACGGCACCCCGGTCACATTACGACAAACCCTGTCACCGGTACG ++ +ADG69ADG69ADG69ADG69ADG69 +@one_mate_reads_through +TCCGGCCCCCCAGATGGTTGTAGCAATACCGGAAA ++ +ADG69 +@terminal_mismatch +GTGCCGGGCTCAGAGAAGCCGCTTCATGGTTCGACAAGGCTCCTCCCTCAATTAAGATTC ++ +ADG69ADG69ADG69 +@over_max_insert +AATTATCCTTCAAACGTAGGCTGCCGGGCTCAGAGAAGCCGCTTCATGGT ++ +ADG69ADG69ADG69ADG69ADG69ADG69ADG69 +@plain +CACGGATGCAGTACAAAGTTGTATTCTATCAGTAGAGATATTCAACGTACCTCCCGTACG ++ +ADG69ADG69ADG69 diff --git a/crates/eschalign-cli/tests/fixtures/pairs/pairs_panel.fa b/crates/eschalign-cli/tests/fixtures/pairs/pairs_panel.fa new file mode 100644 index 0000000..41d7861 --- /dev/null +++ b/crates/eschalign-cli/tests/fixtures/pairs/pairs_panel.fa @@ -0,0 +1,8 @@ +>P1 +CAGTGCCTCTGCCACCGTGCCCACTTCGACGTGGGGAAGCCGTACCGGTGACAGGGTTTGTCGTAATGTGACCGGGGTGCCGTTTCCATCACACGACAGCAGTTCCGGAGAGACCATAATGGGGCGTTTTGTATATACGCCCCGTAGATT +>P2 +CAGTGCCTCTGCCACCGTGCCCACTTCGACGTGGGGAAGCCGTACCGGTGACAGGGTTTGTCGTAATGTGACCGGGGTGCCGTTTCCATCACACGACAGCCGTTCCGGAGAGACCATAATGGGGCGTTTTGTATATACGCCCCGTAGATT +>U1 +GTCACAAGGAGGTAGCTGTGCCAATTCGTCATGAGTGAGCCATCCAACATACGTGGCCGGCCGGGTCTAGGAATCTTAATTGAGGGAGGAGCCTTGTCGAACCATGAAGCGGCTTCTCTGAGCCCGGCAGCCTACGTTTGAAGGATAATT +>U2 +TTTTGTCGACAAACGTCCGTCTTGCGTATGGCACCTTTCCGGTATTGCTACAACCATCTGGGGGGCCGGATCTATTCTGGCGTACGGGAGGTACGTTGAATATCTCTACTGATAGAATACAACTTTGTACTGCATCCGTGGCAAACTAAC diff --git a/crates/eschalign-cli/tests/fixtures/unpaired_golden.tsv b/crates/eschalign-cli/tests/fixtures/unpaired_golden.tsv new file mode 100644 index 0000000..11720e6 --- /dev/null +++ b/crates/eschalign-cli/tests/fixtures/unpaired_golden.tsv @@ -0,0 +1,6 @@ +default 60 ba7eed30a28846b1 +semiglobal_both_secondary 63 dd74fd88b4017599 +sort_coordinate 60 fe04859988b79c33 +min_score_max_ties 60 ba7eed30a28846b1 +scoring_bwa 60 0286fbe4cd147248 +max_read_len 60 5829bf43aa120d38 diff --git a/crates/eschalign/CLAUDE.md b/crates/eschalign/CLAUDE.md index e4cd5d0..eccfec2 100644 --- a/crates/eschalign/CLAUDE.md +++ b/crates/eschalign/CLAUDE.md @@ -11,9 +11,10 @@ - `scoring.rs` — `Scoring`/`Mode`: a gap of length *k* costs `gap_open + (k-1)·gap_extend` (parasail's convention with signs, so bwa's `-O1 -E1` is `-2,-1` here); `Local` = Smith–Waterman, `SemiGlobal` = overlap (leading/trailing gaps of *either* sequence free = parasail `sg`). - `scalar.rs` — the Gotoh oracle: score-only and with traceback. Every other path is tested against it by equality — integer scoring throughout, no tolerance anywhere. - `simd/` — `avx2.rs`/`avx512.rs` score-only inter-sequence kernels (one reference per i16 lane) in two loop orders — `score_group` (row-major) and `score_group_transposed` (reference-major, #415) — plus `trace_pairs` (one (read, reference) pair per lane); `mod.rs` owns `Backend` (`available`, `best`, `cap_from_env`), `Kernel`, `TRANSPOSED_MIN_READ_LEN`, `fits_i16` and the `Profile` both score kernels read (`prof` row-major, `tprof` code-major with −∞ padding, built lazily from `prof` in a `OnceLock` the first time the transposed kernel runs on a group, `ends` per column). +- `adapter.rs` — 3' adapter read-through (`find`, `detect_pair`): a pure function the CLI calls before aligning a read pair; the adapter is *not* aligned, so terminal mismatches stay mismatches. The mates must agree where the insert ends (`len - ad` equal) or the weaker detection is dropped. - `sam.rs` — CIGAR strings and the `MD`/`NM` pair, reproducing `samtools calmd` including at reference ambiguity codes. - `panel.rs` — `Panel`: the reference set in file order, duplicate names refused. -- `mapper.rs` — `Aligner` (panel + dispatch + tie rules): `map_reads` (CPU, scores then traces) and `map_reads_scored` (takes a precomputed `ScoreMatrix`, e.g. from the GPU, and does everything after scoring identically); `kernel_for` (which loop order a read length gets; `with_transposed_min_len` moves or disables the switch — the A/B lever); `score_groups`/`score_group` (one lane group at a time, so a caller can spread one long read over its threads — the crate itself spawns nothing); `Hit`, `ReadMapping`, `MapOptions`. +- `mapper.rs` — `Aligner` (panel + dispatch + tie rules): `map_reads` (CPU, scores then traces) and `map_reads_scored` (takes a precomputed `ScoreMatrix`, e.g. from the GPU, and does everything after scoring identically); `kernel_for` (which loop order a read length gets; `with_transposed_min_len` moves or disables the switch — the A/B lever); `score_groups`/`score_group` (one lane group at a time, so a caller can spread one long read over its threads — the crate itself spawns nothing); `Hit`, `ReadMapping`, `MapOptions`. `mapper/pair.rs` — `map_pairs`/`map_pairs_scored`, `PairOptions`, `PairMapping::{Proper, Separate}`, `PairHit`: candidate = (reference, template strand), score = both mates' scores summed with both ≥ `min_score`, traced a score tier at a time (all open pairs' next tier in one `align_pairs` batch) until a tier passes the orientation/insert constraints on the traced coordinates; the tie set is that tier. No candidate → each mate as `map_reads_scored` with both strands. - `cuda/` (`gpu`) — `GpuScorer`: the same score, one CUDA thread per (read, reference) pair, batched. `mod.rs` orchestration, `kernel.rs` the CUDA source string. - `error.rs` — `AlignError` (`FromStr` for `Scoring`/`Mode` parse errors, surfaced by the CLI's `value_parser`). diff --git a/crates/eschalign/src/adapter.rs b/crates/eschalign/src/adapter.rs new file mode 100644 index 0000000..ce1610a --- /dev/null +++ b/crates/eschalign/src/adapter.rs @@ -0,0 +1,198 @@ +// SPDX-License-Identifier: MIT + +//! 3' adapter read-through, found before alignment so that nothing is clipped +//! by the aligner. +//! +//! A mate that reads past the end of a short insert carries the sequencing +//! adapter on its 3' end. Left in, the adapter either costs score against the +//! reference or, in `local` mode, is soft-clipped together with whatever +//! mismatches sit beside it — which would erase the misincorporation and +//! RT-stop signal at the read end. So the adapter is found *first*, by its +//! known prefix, and the caller aligns only the bases before it (end to end on +//! the reference) and writes the adapter bases as a soft clip it knows the +//! length of. +//! +//! # The match rule +//! +//! The adapter's prefix is searched at every position `p` of the read, leftmost +//! first. The compared stretch is `read[p..]` against the adapter's first +//! `min(adapter length, read length - p)` bases, so a read that runs out +//! mid-adapter still counts. A stretch is the adapter when it is at least +//! `min_overlap` long and has at most one mismatch per ten bases of overlap +//! (cutadapt's default error rate, `overlap / 10`, rounded down — so none below +//! ten bases). Letters outside A/C/G/T are wildcards, as everywhere in this +//! crate. Everything from `p` to the end of the read is adapter. +//! +//! # Mates agree where the insert ends +//! +//! Two mates of one molecule that both read through must cut the insert at the +//! same place: `len1 - ad1 == len2 - ad2`. [`detect_pair`] checks it, and where +//! the two disagree one detection is chance. The one with the shorter adapter +//! overlap (the weaker evidence) is dropped; on a tie both are. The check is +//! made in insert coordinates — the read length less the adapter — rather +//! than after alignment: a trimmed mate's 3' end *is* the insert end, so the +//! two are the same quantity, and checking before alignment keeps the +//! alignment's input fixed. + +use crate::alphabet::{WILDCARD, encode}; + +/// TruSeq Read 1 adapter, as it appears on the 3' end of a read-through mate 1. +pub const TRUSEQ_R1: &[u8] = b"AGATCGGAAGAGCACACGTCTGAACTCCAGTCA"; +/// TruSeq Read 2 adapter, as it appears on the 3' end of a read-through mate 2. +pub const TRUSEQ_R2: &[u8] = b"AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT"; +/// Default shortest stretch counted as adapter. +pub const DEFAULT_MIN_OVERLAP: usize = 5; + +/// What to look for. +#[derive(Debug, Clone, Copy)] +pub struct AdapterConfig<'a> { + /// Adapter prefix on mate 1's 3' end. + pub r1: &'a [u8], + /// Adapter prefix on mate 2's 3' end. + pub r2: &'a [u8], + /// Shortest match counted as adapter (bases). + pub min_overlap: usize, +} + +impl Default for AdapterConfig<'static> { + fn default() -> Self { + Self { + r1: TRUSEQ_R1, + r2: TRUSEQ_R2, + min_overlap: DEFAULT_MIN_OVERLAP, + } + } +} + +/// How many bases at the 3' end of `read` are `adapter` (0 = none). +pub fn find(read: &[u8], adapter: &[u8], min_overlap: usize) -> usize { + let read = encode(read); + let adapter = encode(adapter); + let min_overlap = min_overlap.max(1); + if read.len() < min_overlap || adapter.is_empty() { + return 0; + } + for p in 0..=read.len() - min_overlap { + let overlap = adapter.len().min(read.len() - p); + if overlap < min_overlap { + break; + } + let mismatches = (0..overlap) + .filter(|&k| { + let (a, b) = (read[p + k], adapter[k]); + a != b && a != WILDCARD && b != WILDCARD + }) + .count(); + if mismatches <= overlap / 10 { + return read.len() - p; + } + } + 0 +} + +/// Adapter bases at the 3' end of each mate, with the two mates' insert ends +/// reconciled (see the module docs): `[ad1, ad2]`. +pub fn detect_pair(r1: &[u8], r2: &[u8], cfg: &AdapterConfig<'_>) -> [usize; 2] { + let mut ad = [ + find(r1, cfg.r1, cfg.min_overlap), + find(r2, cfg.r2, cfg.min_overlap), + ]; + let (l1, l2) = (r1.len(), r2.len()); + if ad[0] > 0 && ad[1] > 0 && l1 - ad[0] != l2 - ad[1] { + let (e1, e2) = (ad[0].min(cfg.r1.len()), ad[1].min(cfg.r2.len())); + match e1.cmp(&e2) { + std::cmp::Ordering::Greater => ad[1] = 0, + std::cmp::Ordering::Less => ad[0] = 0, + std::cmp::Ordering::Equal => ad = [0, 0], + } + } + ad +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn full_and_partial_adapters_are_found() { + let mut r = b"ACGTACGTAC".to_vec(); + r.extend_from_slice(&TRUSEQ_R1[..20]); + assert_eq!(find(&r, TRUSEQ_R1, 5), 20); + // The read ends inside the adapter's first bases. + let mut r = b"ACGTACGTAC".to_vec(); + r.extend_from_slice(&TRUSEQ_R1[..6]); + assert_eq!(find(&r, TRUSEQ_R1, 5), 6); + // A full adapter followed by more sequence: all of it is adapter. + let mut r = b"ACGTACGTAC".to_vec(); + r.extend_from_slice(TRUSEQ_R1); + r.extend_from_slice(b"TTTT"); + assert_eq!(find(&r, TRUSEQ_R1, 5), TRUSEQ_R1.len() + 4); + } + + #[test] + fn mismatch_tolerance_scales_with_overlap() { + let mut r = b"GGGGGGGGGGGG".to_vec(); + let mut a = TRUSEQ_R1[..20].to_vec(); + a[5] = b'T'; // 1 mismatch in 20 bases: allowed (2) + a[12] = b'T'; // 2 mismatches: still allowed + r.extend_from_slice(&a); + assert_eq!(find(&r, TRUSEQ_R1, 5), 20); + a[15] = b'C'; // 3 in 20: refused + let mut r = b"GGGGGGGGGGGG".to_vec(); + r.extend_from_slice(&a); + assert_eq!(find(&r, TRUSEQ_R1, 5), 0); + // Under ten bases of overlap there is no tolerance at all. + let mut r = b"GGGGGGGGGGGG".to_vec(); + r.extend_from_slice(b"AGATCGGAT"); + assert_eq!(find(&r, TRUSEQ_R1, 5), 0); + } + + #[test] + fn short_overlap_below_the_minimum_is_not_adapter() { + let mut r = b"GGGGGGGGGGGG".to_vec(); + r.extend_from_slice(b"AGAT"); + assert_eq!(find(&r, TRUSEQ_R1, 5), 0); + assert_eq!(find(&r, TRUSEQ_R1, 4), 4); + assert_eq!(find(b"AG", TRUSEQ_R1, 5), 0); + assert_eq!(find(b"", TRUSEQ_R1, 5), 0); + } + + #[test] + fn wildcards_match() { + let mut r = b"GGGGGGGGGGGG".to_vec(); + r.extend_from_slice(b"AGATNGGAAGAG"); + assert_eq!(find(&r, TRUSEQ_R1, 5), 12); + } + + #[test] + fn mates_must_agree_where_the_insert_ends() { + let cfg = AdapterConfig::default(); + let insert = b"GGGGCCCCGGGGCCCCGGGG"; // 20 nt + let mk = |ad: &[u8], extra: usize| { + let mut r = insert.to_vec(); + r.extend_from_slice(&ad[..extra]); + r + }; + // Both read through by 12 bases: they agree. + let (a, b) = (mk(TRUSEQ_R1, 12), mk(TRUSEQ_R2, 12)); + assert_eq!(detect_pair(&a, &b, &cfg), [12, 12]); + // Mate 2's insert ends two bases later: the weaker detection goes. + let shifted = |ad: &[u8], extra: usize| { + let mut r = insert.to_vec(); + r.extend_from_slice(b"CC"); + r.extend_from_slice(&ad[..extra]); + r + }; + let b2 = shifted(TRUSEQ_R2, 9); + assert_eq!(detect_pair(&a, &b2, &cfg), [12, 0]); + let a2 = shifted(TRUSEQ_R1, 7); + assert_eq!(detect_pair(&a2, &b, &cfg), [0, 12]); + // Equal evidence, different ends: neither. + let a3 = mk(TRUSEQ_R1, 8); + let b3 = shifted(TRUSEQ_R2, 8); + assert_eq!(detect_pair(&a3, &b3, &cfg), [0, 0]); + // One mate only: stands. + let clean = insert.to_vec(); + assert_eq!(detect_pair(&a, &clean, &cfg), [12, 0]); + } +} diff --git a/crates/eschalign/src/lib.rs b/crates/eschalign/src/lib.rs index 7703bb4..17595d0 100644 --- a/crates/eschalign/src/lib.rs +++ b/crates/eschalign/src/lib.rs @@ -30,6 +30,8 @@ //! i16 saturating), AVX2 and AVX-512BW, runtime-dispatched. //! * [`sam`] — CIGAR strings and the `MD`/`NM` pair, reproducing //! `samtools calmd` including at reference ambiguity codes. +//! * [`adapter`] — 3' adapter read-through detection, run before a read pair +//! is aligned. //! * [`Aligner`] — the panel plus the dispatch plus winner selection: what a //! caller actually uses. //! * [`fasta`] — the one reference FASTA parser (rnabioco/escapepod-rs#410): @@ -46,6 +48,7 @@ //! gzipped) FASTA file before handing it to [`fasta::read_fasta`], is //! `eschalign-cli`'s job. +pub mod adapter; pub mod alphabet; #[cfg(feature = "gpu")] pub mod cuda; @@ -59,7 +62,9 @@ mod scoring; pub mod simd; pub use error::AlignError; -pub use mapper::{Aligner, Hit, MapOptions, ReadMapping, ScoreMatrix}; +pub use mapper::{ + Aligner, Hit, MapOptions, PairHit, PairMapping, PairOptions, ReadMapping, ScoreMatrix, +}; pub use panel::{Panel, Reference}; pub use scalar::{Alignment, CigarOp}; pub use scoring::{Mode, Scoring}; diff --git a/crates/eschalign/src/mapper.rs b/crates/eschalign/src/mapper.rs index f51f38a..615edc7 100644 --- a/crates/eschalign/src/mapper.rs +++ b/crates/eschalign/src/mapper.rs @@ -23,6 +23,9 @@ use crate::scalar::Alignment; use crate::scoring::{Mode, Scoring}; use crate::simd::{self, Backend, Kernel, Profile}; +mod pair; +pub use pair::{PairHit, PairMapping, PairOptions}; + /// Per-read policy. The default: minimum score 0, every tie traced, forward /// strand only. #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] diff --git a/crates/eschalign/src/mapper/pair.rs b/crates/eschalign/src/mapper/pair.rs new file mode 100644 index 0000000..a218757 --- /dev/null +++ b/crates/eschalign/src/mapper/pair.rs @@ -0,0 +1,574 @@ +// SPDX-License-Identifier: MIT + +//! Pair-aware winner selection: [`Aligner::map_pairs`]. +//! +//! # Semantics +//! +//! A *candidate* is a reference together with a template strand. On the +//! forward template strand mate 1 is read as sequenced and mate 2 is +//! reverse-complemented (an `FR` pair, mate 1 upstream); with +//! [`PairOptions::both_strands`] the reverse template strand is a candidate +//! too (mate 1 reverse-complemented, mate 2 as sequenced). A candidate's score +//! is `score(mate 1) + score(mate 2)`, each against the reference, and it is +//! eligible only when **both** mates are at or above +//! [`PairOptions::min_score`] — the sum alone does not let one good mate carry +//! a bad one. +//! +//! Eligible candidates are taken in descending score. The first score for +//! which at least one candidate also satisfies the pair constraints is the +//! pair's score, and **all** such candidates at that score are the *tie set* +//! (reference order, forward template strand first); its first member is the +//! primary. The constraints are checked on the traced alignments, so a score +//! that only exists as an impossible pair (mates 10 kb apart, facing away from +//! each other) falls through to the next score rather than winning: +//! +//! * both mates align at least one base; +//! * the forward-strand mate starts before the reverse-strand mate ends (no +//! "outie" — a pair facing away from each other); +//! * the insert — from the leftmost aligned base of either mate to the +//! rightmost — is at most [`PairOptions::max_insert`]. An insert shorter than +//! the reads (the mates overlap completely, one reading through into +//! adapter) is valid, and is what it is: the span, not the read length. +//! +//! If no candidate qualifies the pair is [`PairMapping::Separate`]: each mate +//! is mapped exactly as [`Aligner::map_reads`] maps a read, on **both** +//! strands (a mate that is the wrong way round must be reported where it is, +//! not forced onto the strand the pair rule wanted) with its own tie set. +//! +//! `XS` per mate is that mate's best score over the candidates *outside* the +//! pair's tie set, mirroring the unpaired definition. + +use super::*; + +/// Per-pair policy. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct PairOptions { + /// Both mates must score at least this to be a candidate. + pub min_score: i32, + /// How many tie-set members *besides the primary* are reported + /// (`None` = all). The tie set is always counted in full. + pub max_ties: Option, + /// Longest allowed insert, bases. + pub max_insert: usize, + /// Also consider the reverse template strand (mate 1 reversed). + pub both_strands: bool, +} + +impl Default for PairOptions { + fn default() -> Self { + Self { + min_score: 0, + max_ties: None, + max_insert: 500, + both_strands: false, + } + } +} + +/// One traced placement of a pair on one reference. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct PairHit { + /// Panel index of the reference. + pub reference: usize, + /// The template is on the reverse strand: mate 1 aligned reverse- + /// complemented, mate 2 as sequenced. + pub reverse: bool, + /// Mate 1's and mate 2's alignments. Each [`Hit::reverse`] is that mate's + /// own strand, and its coordinates are on the mate as aligned. + pub mates: [Hit; 2], + /// Leftmost aligned reference base of either mate to the rightmost. + pub insert: usize, +} + +impl PairHit { + /// The pair score: the two mates' scores summed. + pub fn score(&self) -> i32 { + self.mates[0].alignment.score + self.mates[1].alignment.score + } +} + +/// What one read pair produced. +#[derive(Debug, Clone, PartialEq, Eq)] +pub enum PairMapping { + /// A pair that satisfies every constraint on some reference. + Proper { + /// The pair score (both mates' scores summed). + score: i32, + /// Size of the tie set (≥ 1). + n_tied: usize, + /// Per mate: the best score outside the tie set (SAM `XS`). + suboptimal: [Option; 2], + /// Primary first, then up to `max_ties` more tie-set members. + hits: Vec, + }, + /// No reference satisfies the constraints: each mate at its own best + /// reference, as an unpaired read, on either strand. + Separate([ReadMapping; 2]), +} + +/// One eligible candidate. +#[derive(Clone, Copy)] +struct Cand { + sum: i32, + reference: usize, + reverse: bool, + scores: [i32; 2], +} + +struct PairState<'a> { + letters: [&'a [u8]; 2], + /// `[mate][reverse]`: codes of the mate as sequenced / reverse-complemented. + codes: [[Vec; 2]; 2], + /// `[mate][reverse]`: the panel scores, when this pair has candidates. + rows: [[Vec; 2]; 2], + cands: Vec, + /// `cands` ranges of equal sum, best first. + tiers: Vec<(usize, usize)>, + next_tier: usize, + /// The tie set, once found. + found: Vec<(Cand, [Alignment; 2])>, +} + +impl PairState<'_> { + /// Mate `m`'s orientation in a candidate on template strand `reverse`: + /// mate 1 follows the template, mate 2 is its opposite. + fn mate_reverse(m: usize, template_reverse: bool) -> bool { + template_reverse ^ (m == 1) + } + + fn open(&self) -> bool { + self.found.is_empty() && self.next_tier < self.tiers.len() + } +} + +impl Aligner { + /// Align read pairs (letters, as sequenced; mate 1 then mate 2), one + /// [`PairMapping`] each, in order. See the module docs for the rules. + pub fn map_pairs(&self, pairs: &[[&[u8]; 2]], opts: &PairOptions) -> Vec { + self.map_pairs_scored(pairs, opts, |_, _, _| None) + } + + /// [`Aligner::map_pairs`] with some panel scores already computed + /// elsewhere (the CUDA kernel). `prescored(pair, mate, reverse)` is that + /// mate's score against every reference, in panel order, as sequenced + /// (`reverse == false`) or reverse-complemented; `None` scores it here. + /// The scores are the caller's contract, exactly as for + /// [`Aligner::map_reads_scored`]: they must equal [`Aligner::score_all`]'s. + pub fn map_pairs_scored<'s>( + &self, + pairs: &[[&[u8]; 2]], + opts: &PairOptions, + prescored: impl Fn(usize, usize, bool) -> Option<&'s [i16]>, + ) -> Vec { + let n_refs = self.panel.len(); + let strands: &[bool] = if opts.both_strands { + &[false, true] + } else { + &[false] + }; + + // Phase 1: score, and rank the eligible candidates. + let mut states: Vec> = pairs + .iter() + .enumerate() + .map(|(p, letters)| { + let fwd = [encode(letters[0]), encode(letters[1])]; + let rev = [ + alphabet::reverse_complement_codes(&fwd[0]), + alphabet::reverse_complement_codes(&fwd[1]), + ]; + let mut st = PairState { + letters: *letters, + codes: [ + [fwd[0].clone(), rev[0].clone()], + [fwd[1].clone(), rev[1].clone()], + ], + rows: Default::default(), + cands: Vec::new(), + tiers: Vec::new(), + next_tier: 0, + found: Vec::new(), + }; + if fwd[0].is_empty() || fwd[1].is_empty() { + return st; + } + for &t in strands { + for m in 0..2 { + let rv = PairState::mate_reverse(m, t); + if st.rows[m][usize::from(rv)].is_empty() { + let mut row = Vec::new(); + self.score_or_take( + &st.codes[m][usize::from(rv)], + prescored(p, m, rv), + &mut row, + ); + st.rows[m][usize::from(rv)] = row; + } + } + let r0 = &st.rows[0][usize::from(PairState::mate_reverse(0, t))]; + let r1 = &st.rows[1][usize::from(PairState::mate_reverse(1, t))]; + for r in 0..n_refs { + let (a, b) = (r0[r], r1[r]); + if a >= opts.min_score && b >= opts.min_score { + st.cands.push(Cand { + sum: a + b, + reference: r, + reverse: t, + scores: [a, b], + }); + } + } + } + // Best sum first; ties in reference order, forward strand first. + st.cands + .sort_by_key(|c| (std::cmp::Reverse(c.sum), c.reference, c.reverse)); + let mut start = 0; + for k in 1..=st.cands.len() { + if k == st.cands.len() || st.cands[k].sum != st.cands[start].sum { + st.tiers.push((start, k)); + start = k; + } + } + st + }) + .collect(); + + // Phase 2: trace tier by tier, every open pair's next tier together, + // until each pair has a tie set or runs out of tiers. + loop { + let plan: Vec<(usize, usize, usize)> = states + .iter() + .enumerate() + .filter(|(_, st)| st.open()) + .map(|(p, st)| { + let (a, b) = st.tiers[st.next_tier]; + (p, a, b) + }) + .collect(); + if plan.is_empty() { + break; + } + let traced = { + let mut jobs: Vec<(&[u8], &[u8])> = Vec::new(); + for &(p, a, b) in &plan { + let st = &states[p]; + for c in &st.cands[a..b] { + let reference = &self.panel.get(c.reference).codes[..]; + for m in 0..2 { + let rv = PairState::mate_reverse(m, c.reverse); + jobs.push((&st.codes[m][usize::from(rv)], reference)); + } + } + } + simd::align_pairs(self.backend, &jobs, &self.scoring, self.mode) + }; + let mut traced = traced.into_iter(); + for &(p, a, b) in &plan { + let st = &mut states[p]; + let mut found = Vec::new(); + for k in a..b { + let c = st.cands[k]; + let al = [traced.next().unwrap(), traced.next().unwrap()]; + debug_assert!( + al[0].score == c.scores[0] && al[1].score == c.scores[1], + "traceback disagrees with the kernel" + ); + if self.pair_ok(&al, c.reverse, opts) { + found.push((c, al)); + } + } + if found.is_empty() { + st.next_tier += 1; + } else { + st.found = found; + } + } + } + + // Phase 3: build the proper pairs; collect the rest for the fallback. + let mut out: Vec> = Vec::with_capacity(states.len()); + let mut separate: Vec = Vec::new(); + for (p, st) in states.iter_mut().enumerate() { + if st.found.is_empty() { + separate.push(p); + out.push(None); + continue; + } + let found = std::mem::take(&mut st.found); + let n_tied = found.len(); + let score = found[0].0.sum; + // Best score outside the tie set, per mate. + let mut in_set = vec![false; 2 * n_refs]; + for (c, _) in &found { + in_set[2 * c.reference + usize::from(c.reverse)] = true; + } + let mut suboptimal = [None; 2]; + for &t in strands { + for r in 0..n_refs { + if in_set[2 * r + usize::from(t)] { + continue; + } + for (m, sub) in suboptimal.iter_mut().enumerate() { + let s = st.rows[m][usize::from(PairState::mate_reverse(m, t))][r]; + *sub = Some(sub.map_or(s, |x: i32| x.max(s))); + } + } + } + let keep = opts + .max_ties + .map_or(n_tied, |k| n_tied.min(k.saturating_add(1))); + let rc = [ + alphabet::reverse_complement(st.letters[0]), + alphabet::reverse_complement(st.letters[1]), + ]; + let hits = found + .into_iter() + .take(keep) + .map(|(c, al)| { + let reference = self.panel.get(c.reference); + let [a0, a1] = al; + let mut mates = [a0, a1].into_iter().enumerate().map(|(m, alignment)| { + let rv = PairState::mate_reverse(m, c.reverse); + let letters: &[u8] = if rv { &rc[m] } else { st.letters[m] }; + let (md, nm) = sam::md_nm( + letters, + &reference.seq, + alignment.query_start, + alignment.ref_start, + &alignment.ops, + ); + Hit { + reference: c.reference, + reverse: rv, + alignment, + md, + nm, + } + }); + let mates = [mates.next().unwrap(), mates.next().unwrap()]; + let insert = span(&mates); + PairHit { + reference: c.reference, + reverse: c.reverse, + mates, + insert, + } + }) + .collect(); + out.push(Some(PairMapping::Proper { + score, + n_tied, + suboptimal, + hits, + })); + } + + // Pairs with no qualifying candidate: each mate on its own, either strand. + if !separate.is_empty() { + let flat: Vec<&[u8]> = separate + .iter() + .flat_map(|&p| [pairs[p][0], pairs[p][1]]) + .collect(); + let each = MapOptions { + min_score: opts.min_score, + max_ties: opts.max_ties, + both_strands: true, + }; + let mut mapped = self + .map_reads_scored(&flat, &each, |k, rev| { + prescored(separate[k / 2], k % 2, rev) + }) + .into_iter(); + for &p in &separate { + let a = mapped.next().unwrap(); + let b = mapped.next().unwrap(); + out[p] = Some(PairMapping::Separate([a, b])); + } + } + out.into_iter().map(Option::unwrap).collect() + } + + /// The pair constraints on two traced mates (see the module docs). + fn pair_ok(&self, al: &[Alignment; 2], template_reverse: bool, opts: &PairOptions) -> bool { + if al[0].is_empty() || al[1].is_empty() { + return false; + } + // The mate read as sequenced on the template strand is the "forward" + // one: mate 1 on the forward template strand, mate 2 on the reverse. + let (fwd, rev) = if template_reverse { + (&al[1], &al[0]) + } else { + (&al[0], &al[1]) + }; + if fwd.ref_start >= rev.ref_end { + return false; + } + let lo = al[0].ref_start.min(al[1].ref_start); + let hi = al[0].ref_end.max(al[1].ref_end); + hi - lo <= opts.max_insert + } +} + +/// Leftmost aligned reference base to rightmost, over both mates. +fn span(mates: &[Hit; 2]) -> usize { + let lo = mates[0] + .alignment + .ref_start + .min(mates[1].alignment.ref_start); + let hi = mates[0].alignment.ref_end.max(mates[1].alignment.ref_end); + hi - lo +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Two paralogs that differ at one position, 60 nt each, in a panel with + /// an unrelated third reference. + fn panel() -> (Panel, Vec, Vec) { + let mut x = 0x2545_f491_4f6c_dd1du64; + let mut base = || { + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + b"ACGT"[(x % 4) as usize] + }; + let a: Vec = (0..60).map(|_| base()).collect(); + let mut b = a.clone(); + b[50] = if a[50] == b'A' { b'C' } else { b'A' }; + let c: Vec = (0..60).map(|_| base()).collect(); + let p = Panel::new([("a", a.clone()), ("b", b.clone()), ("c", c)]).unwrap(); + (p, a, b) + } + + fn aligner(mode: Mode) -> (Aligner, Vec, Vec) { + let (p, a, b) = panel(); + (Aligner::new(p, Scoring::default(), mode).unwrap(), a, b) + } + + #[test] + fn pair_picks_the_paralog_its_mates_support() { + let (al, a, _) = aligner(Mode::SemiGlobal); + // Mate 1 covers 0..30 (identical in both paralogs); mate 2 covers + // 30..60 of the template, reverse-complemented, and carries a's base. + let r1 = a[..30].to_vec(); + let r2 = alphabet::reverse_complement(&a[30..]); + let alone = al.map_read(&r1, &MapOptions::default()); + assert!(alone.n_tied >= 2, "mate 1 alone is a tie"); + let m = al.map_pairs(&[[&r1, &r2]], &PairOptions::default()); + match &m[0] { + PairMapping::Proper { + n_tied, + hits, + score, + .. + } => { + assert_eq!(*n_tied, 1); + assert_eq!(hits[0].reference, 0); + assert_eq!(*score, 120); + assert_eq!(hits[0].insert, 60); + assert!(!hits[0].mates[0].reverse && hits[0].mates[1].reverse); + } + other => panic!("not proper: {other:?}"), + } + } + + #[test] + fn over_max_insert_and_wrong_orientation_are_separate() { + let (al, a, _) = aligner(Mode::SemiGlobal); + let r1 = a[..20].to_vec(); + let r2 = alphabet::reverse_complement(&a[40..]); + let opts = PairOptions { + max_insert: 50, + min_score: 30, + ..PairOptions::default() + }; + assert!(matches!( + al.map_pairs(&[[&r1, &r2]], &opts)[0], + PairMapping::Separate(_) + )); + // Semi-global scoring gives any read some positive score somewhere, + // so a floor is what makes a wrong-way mate "not aligned". + let roomy = PairOptions { + min_score: 30, + ..PairOptions::default() + }; + assert!(matches!( + al.map_pairs(&[[&r1, &r2]], &roomy)[0], + PairMapping::Proper { .. } + )); + // Both mates as sequenced on the same strand: FF. + let ff2 = a[40..].to_vec(); + assert!(matches!( + al.map_pairs(&[[&r1, &ff2]], &roomy)[0], + PairMapping::Separate(_) + )); + // An outie: mate 2 entirely upstream of mate 1. + let up = alphabet::reverse_complement(&a[..20]); + let down = a[40..].to_vec(); + assert!(matches!( + al.map_pairs(&[[&down, &up]], &roomy)[0], + PairMapping::Separate(_) + )); + } + + #[test] + fn reverse_template_strand_only_with_both_strands() { + let (al, a, _) = aligner(Mode::SemiGlobal); + // Mate 1 reverse-complemented from the downstream part, mate 2 forward upstream. + let r1 = alphabet::reverse_complement(&a[40..]); + let r2 = a[..20].to_vec(); + let fwd = PairOptions { + min_score: 30, + ..PairOptions::default() + }; + assert!(matches!( + al.map_pairs(&[[&r1, &r2]], &fwd)[0], + PairMapping::Separate(_) + )); + let both = PairOptions { + both_strands: true, + ..fwd + }; + match &al.map_pairs(&[[&r1, &r2]], &both)[0] { + PairMapping::Proper { hits, .. } => { + assert!(hits[0].reverse); + assert!(hits[0].mates[0].reverse && !hits[0].mates[1].reverse); + } + other => panic!("not proper: {other:?}"), + } + } + + #[test] + fn both_mates_must_clear_min_score() { + let (al, a, _) = aligner(Mode::SemiGlobal); + let r1 = a[..30].to_vec(); + let r2 = alphabet::reverse_complement(&a[30..]); + let strict = PairOptions { + min_score: 61, + ..PairOptions::default() + }; + // Each mate scores 60: the sum (120) clears 61, a mate (60) does not. + assert!(matches!( + al.map_pairs(&[[&r1, &r2]], &strict)[0], + PairMapping::Separate(_) + )); + } + + #[test] + fn batched_pairs_equal_one_at_a_time() { + let (al, a, b) = aligner(Mode::SemiGlobal); + let reads: Vec<[Vec; 2]> = vec![ + [a[..30].to_vec(), alphabet::reverse_complement(&a[30..])], + [b[..30].to_vec(), alphabet::reverse_complement(&b[30..])], + [a[..20].to_vec(), a[40..].to_vec()], + [Vec::new(), a[..20].to_vec()], + ]; + let pairs: Vec<[&[u8]; 2]> = reads.iter().map(|r| [&r[0][..], &r[1][..]]).collect(); + let opts = PairOptions::default(); + let all = al.map_pairs(&pairs, &opts); + for (p, m) in pairs.iter().zip(&all) { + assert_eq!(&al.map_pairs(&[*p], &opts)[0], m); + } + } +} diff --git a/docs/align.md b/docs/align.md index 3d9bdc2..35a915b 100644 --- a/docs/align.md +++ b/docs/align.md @@ -29,7 +29,7 @@ eschalign [OPTIONS] -r -o | Argument | Description | |----------|-------------| -| `` | Reads: an unaligned BAM (dorado's uBAM), or FASTA/FASTQ, plain or gzip. The format is sniffed from the content, not the extension | +| `` | Reads: an unaligned BAM (dorado's uBAM), or FASTA/FASTQ, plain or gzip. The format is sniffed from the content, not the extension. Omitted with `-1`/`-2` | A BAM whose records are already aligned is accepted too, as long as they are on the forward strand (the read as sequenced): the alignment fields and any @@ -46,7 +46,7 @@ reverse-strand record is refused. | `--sort ` | Record order (default `unsorted`) — see [Sorted output](#sorted-output) | | `--sort-memory ` | With `--sort coordinate`, records held in memory before spilling to a temporary file (default `4G`; `K`/`M`/`G`/`T`, binary) | | `--tmp-dir ` | With `--sort coordinate`, where the spill file goes (default: the output's directory; `$TMPDIR` for `-o -`) | -| `--mode ` | `local` (default) or `semiglobal` — see [Modes](#modes) | +| `--mode ` | `local` (default; `semiglobal` for paired input) or `semiglobal` — see [Modes](#modes) | | `--scoring ` | Match, mismatch, gap open, gap extend (default `2,-1,-10,-1`) — see [Scoring](#scoring) | | `--min-score ` | A read whose best score is below this is written unmapped (default `0`) | | `--max-read-len ` | A read longer than this (nt) is not aligned at all and is written unmapped (default `1000`; `0` = no limit) — see [Long reads](#long-reads) | @@ -56,6 +56,12 @@ reverse-strand record is refused. | `--read-ids ` | Align only the reads named in FILE, one per line (like `samtools view -N`); the rest are not written | | `--batch-size ` | Reads per batch (default `50000`); memory is proportional to it | | `-t, --threads ` | Threads for parallel processing | +| `-1, --r1 `, `-2, --r2 ` | Paired input: mate 1 and mate 2 files (plain or gzip). Mates are matched by name, `/1` and `/2` stripped — see [Paired-end](#paired-end) | +| `--interleaved` | `READS` is an interleaved FASTQ/FASTA (mate 1, mate 2, …). A BAM's pairing comes from its flags and needs no option | +| `--max-insert ` | Paired only: longest insert a proper pair may have (default `500`) | +| `--adapter-r1 `, `--adapter-r2 ` | Paired only: 3' adapter of each mate (default TruSeq Read 1 / Read 2) | +| `--min-adapter-overlap ` | Paired only: shortest 3' match counted as adapter (default `5`) | +| `--no-adapter-trim` | Paired only: no adapter detection, no `ad` tag | | `--device ` | Where panel scoring runs (a `gpu` build only; tracebacks and output are always on the CPU, and the output is identical either way). `auto` (default) scores on a visible CUDA GPU, else on the CPU; `cpu` forces the CPU; `gpu` requires the GPU and fails rather than falling back. See [GPU scoring](#gpu-scoring) | ## Scoring @@ -180,6 +186,175 @@ one `N` has `NM:i:1` and an `MD` with that `N` in it, and `AS` counts the `N` as a match. A pipeline stage that used `samtools calmd` after bwa sees the same `MD`/`NM` for the same CIGAR. +## Paired-end + +`-1 R1 -2 R2`, `--interleaved`, or a BAM whose first record carries 0x1 with +0x40 or 0x80 switches on paired-end alignment (rnabioco/eschalign#1). The three +input forms give identical output for the same pairs, and the unpaired path — +every run without paired input — is unchanged, byte for byte. + +```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 +``` + +Built for snRNAs up to ~250 nt captured by two overlapping mates, and for any +panel with paralogs: the partner mate resolves a mate that is ambiguous alone, +and the signal at the read ends — RT stops, misincorporation — is kept. Better +handling comes first; there is no speed work in it (a paired run scores every +mate once, and traces a tier of candidates at a time). + +**Use `--min-score`.** Semi-global scoring (the paired default) gives every read +*some* positive score somewhere, so without a floor a mate that matches nothing +is still "aligned" and a wrong-way pair can find a junk partner. The floor +applies to each mate (below). + +### Input + +* **Two files**: `-1 R1.fq -2 R2.fq` (FASTQ or FASTA, plain or gzip); record *k* + of one is the mate of record *k* of the other. +* **Interleaved**: `READS --interleaved`, mate 1 then mate 2. +* **BAM**: an unaligned BAM with 0x1 and 0x40/0x80, mates adjacent. Pairing is + taken from the first primary record, so no option is needed. (A BAM with those + flags used to run as unpaired reads, the flags discarded; it is now a paired + run. A forward-strand aligned BAM is still accepted, as it is for single + reads; its alignment is discarded.) + +Mates are matched on name with a trailing `/1` or `/2` stripped. **An orphan +mate or a name mismatch is a hard error naming the read** — a mate file that is +a record short, an interleaved file with a trailing record, mates that are not +adjacent in a BAM, a paired-flagged record in an otherwise unpaired BAM. Both +mates are written under the stripped name. `--read-ids` selects pairs by that +name. Secondary and supplementary records are skipped as for single reads. + +### Adapters: detected first, bases kept + +`local` mode soft-clips mismatches near a read end, which would erase the +misincorporation and RT-stop signal at the 3' end that this mode is for. So the +adapter is found *before* alignment, on each mate's 3' end: + +1. The adapter's prefix (`--adapter-r1`/`--adapter-r2`; TruSeq Read 1 + `AGATCGGAAGAGCACACGTCTGAACTCCAGTCA` and Read 2 + `AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT`) is searched at every position of the + read, leftmost first. `read[p..]` is compared with the adapter's first + `min(adapter length, bases left)` bases, so a read that ends inside the + adapter still matches. A stretch counts when it is at least + `--min-adapter-overlap` long (default 5) and has at most one mismatch per ten + bases of overlap (`overlap / 10`, rounded down: none below ten bases); `N` and + other non-ACGT letters match anything. +2. The mates must agree where the insert ends. Both reading through means + `len1 − ad1 = len2 − ad2`; where they disagree one detection is chance, so + the one with the shorter adapter overlap is dropped (both, on a tie). The + check is made on the read less its adapter rather than after alignment — a + trimmed mate's 3' end *is* the insert end, so it is the same quantity — which + keeps the alignment's input fixed. +3. Only the bases before the adapter are aligned. The adapter bases are written + as a soft clip (at the end of `SEQ` for a forward mate, the start for a + reverse one), and `ad:i` says how many there were. **Nothing is discarded**: + `SEQ` and `QUAL` are the whole read. + +`ad:i` is each mate's *own* adapter length, and is present (0 when none) on every +primary record of a paired run that searched; with `--no-adapter-trim` it is +absent. When only one mate reads through — a shorter or trimmed partner — that +mate carries the count and the other carries `ad:i:0`. + +### Aligning end to end + +The remainder of each mate is aligned so that a terminal mismatch stays a +mismatch: paired input defaults to `--mode semiglobal`, in which a mate's ends +are kept unless the alignment reaches the end of the reference (where the +overhang is free). A terminal mismatch is then in `MD`/`NM`, not clipped: + +``` +mate 1, last base a mismatch: CIGAR 60M MD:Z:59G0 NM:i:1 (not 59M1S) +``` + +An explicit `--mode local` is honoured — and clips again; a warning says so. + +### Choosing the reference + +A *candidate* is a reference and a template strand. On the forward template +strand mate 1 is read as sequenced and mate 2 is reverse-complemented (an `FR` +pair, mate 1 upstream — the layout of a stranded sense library); `--strand +both` adds the reverse template strand (mate 1 reverse-complemented, mate 2 as +sequenced). Candidate score = `score(mate 1) + score(mate 2)` against the +reference. + +1. A candidate is eligible only if **both** mates are at or above `--min-score`. + (The spec's open question: both, not just the sum — one good mate must not + carry a bad one.) +2. Eligible candidates are taken in descending score, and the first score at + which some candidate also satisfies the constraints wins; **every** candidate + at that score that satisfies them is the tie set (reference order, forward + strand first), its first member the primary. A better-scoring candidate that + violates a constraint falls through to the next score rather than winning. +3. The constraints, on the traced alignments: + * both mates align at least one base; + * orientation: the mate on the template's forward strand starts before the + other mate ends (no "outie"); two mates on the same strand (`FF`/`RR`) are + not a candidate at all; + * insert ≤ `--max-insert` (default 500), the insert being the span from the + leftmost aligned reference base of either mate to the rightmost. An insert + shorter than the read length — full overlap and read-through — is valid. +4. If no candidate qualifies the pair is written **without 0x2**, and each mate + is mapped as an unpaired read would be, on **both** strands — a mate that is + the wrong way round is reported where it is, not forced onto the strand the + pair rule wanted — with its own tie set, MAPQ, `XA` and `XS`. + +### Paired output + +Records are in input order with the two mates adjacent, mate 1 then mate 2; +`--secondary` records follow the pair. With `--sort coordinate` the order is +`samtools sort`'s, so mates are placed by their own coordinates and each record's +`RNEXT`/`PNEXT` still name its partner's place. + +| Field | Paired value | +|-------|--------------| +| `FLAG` | 0x1; 0x2 only for a proper pair; 0x10 / 0x20 for each mate's / its partner's strand; 0x40 / 0x80; 0x4 / 0x8 for an unmapped mate / partner; QC-fail and duplicate (0x200, 0x400) from the input | +| `RNEXT` / `PNEXT` | the partner's reference and position (`=` is the same reference); for an unmapped partner, the mate's own | +| `TLEN` | the true insert — leftmost aligned base to rightmost, over both mates — positive on the leftmost mate (mate 1 if they start together), negative on the other. Set for any pair with both mates on one reference, proper or not | +| `MC:Z` | the partner's CIGAR (absent when the partner is unmapped) | +| `AS:i` | that mate's own score (the pair score is the two summed) | +| `XS:i` | that mate's best score over the candidates outside the pair's tie set | +| `XA:Z` | `ref,±pos,CIGAR,NM;` — see below | +| `NM:i`, `MD:Z` | per mate, as `samtools calmd` | +| `ad:i` | adapter bases at the 3' end — see above | +| MAPQ | **60** for a unique best reference, 0 for a tie — for both mates of a pair, nothing in between | + +An unmapped mate (it matched nothing at `--min-score`) is written the way +`samtools` does: flag 0x4, placed at its mapped partner's reference and +position, and the partner carries 0x8. The paired-run summary reports pairs and +how many were proper. + +**`XA` in a pair.** The open question was whether pair-level `XA` encodes both +mates' positions. It does not: each mate's `XA` lists *that mate's* alignment in +every other tied pair, in bwa's `ref,±pos,CIGAR,NM;` format, and the entries run +in the same order on both mates — entry *k* of mate 1 and entry *k* of mate 2 +are the two halves of tied pair *k*. That keeps `XA` the grammar every parser +already reads, and the partner position of an entry is in the same position of +the other mate's tag. `--secondary` writes the same ties as 0x100 records, mate 1 +then mate 2 for each further pair, whose mate fields point at the partner of +that same pair; no `SEQ`/`QUAL`, tags `NM`, `MD`, `AS`, `RG` only. + +Every input tag is copied byte for byte, as for single reads, on both mates; +`NM`/`MD`/`AS`/`XS`/`XA`/`MC`/`ad` are the aligner's. Mates over `--max-read-len` +follow the existing rule: that mate is written unmapped, and its partner is +aligned alone. + +### GPU + +`--device gpu` scores each mate (as sequenced and reverse-complemented, after +adapter trimming) on the device; pair resolution, tracebacks and records are the +CPU path's own code, so paired output is byte-identical either way. + +### Not in this mode + +Overlap merging into a consensus read, joint paired DP, estimating the insert +distribution from the data, and orphan/single-end handling inside a paired run +(an orphan is an error). The long-read path of `--max-read-len 0` is not used for +mates. + ## Kernels Scoring runs on inter-sequence SIMD kernels — one lane per reference, 16 lanes @@ -308,6 +483,9 @@ eschalign calls.ubam -r trna.fa -o aligned.bam --sort coordinate # Ties as secondary records, at most 5 per read eschalign calls.ubam -r trna.fa -o aligned.bam --secondary --max-ties 5 + +# Illumina read pairs, adapter read-through kept as soft clip (ad:i) +eschalign -1 R1.fq.gz -2 R2.fq.gz -r snrna.fa -o aligned.bam --min-score 40 --max-insert 300 ``` ## Environment variables