eschalign¶

All-against-all alignment of reads to a small reference panel: every read is scored against every reference, the best-scoring reference is reported, and references that tie with it are reported as ties rather than hidden behind a MAPQ of zero. Built for tRNA panels (~150–200 references of ~150 nt), where it replaces bwa mem + samtools calmd.

eschalign reads.ubam -r trna.fa -o aligned.bam

One flat binary, no subcommand. Panel scoring on the GPU needs a build with --features gpu; everything else is in the default build.

There is no seed index. Memory is one batch of reads, and runtime is a function of reads × references × lengths alone — the same for every sample of the same size. That is the right trade for a panel of short references and the wrong one for a genome: the intended range is ≤ ~10 k references of ≤ ~1 kb, and a panel beyond it gets a warning and proportionally more runtime.

Usage¶

eschalign [OPTIONS] -r <FASTA> -o <BAM> <READS>

Arguments¶

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 NM/MD/AS/XS/XA are discarded and the read is aligned afresh. Secondary and supplementary input records are skipped (they are not reads); a reverse-strand record is refused.

Options¶

Scoring¶

--scoring M,X,O,E takes four signed integers. The convention, stated once:

A gap of length k contributes O + (k − 1)·E.

O is therefore the whole cost of a one-base gap. This is parasail's and EMBOSS's convention (with the signs written out), and the aligner is tested against parasail on it. bwa states its penalties the other way round — a gap of length k costs O + k·E — so bwa's -A1 -B1 -O1 -E1 (-x ont2d) is --scoring 1,-1,-2,-1 here.

The default, 2,-1,-10,-1, is gpu-tRNA-mapper's, so a user moving from it sees the same behaviour. It is not a recommendation for any pipeline; choosing a scheme is a validation question (see Replacing bwa).

Every letter other than A/C/G/T/U (upper or lower case; U is T) is a wildcard: pairing it with anything scores as a match. The tRNA panels this is built for carry an N in every reference, and scoring it as a mismatch would penalise every read at the same column for a base the reference does not claim to know.

Modes¶

local — Smith–Waterman. The alignment may start and end anywhere; both overhangs of the read are soft-clipped.

semiglobal — overlap alignment: leading and trailing gaps of either sequence are free. The path starts on the top or left border of the DP matrix and ends on the bottom or right border, and interior gaps are scored normally. In practice: a read that runs off either end of its reference is clipped there for free, and a reference end the read does not reach costs nothing. This is our definition (parasail's sg), not a reproduction of any other tool's semi-global mode.

Optimal alignments are rarely unique. The choice is fixed: the end cell with the smallest reference end, then the smallest read end; from it, the traceback prefers a diagonal step, then an insertion, then a deletion, extends a gap rather than opening one, and (local mode) stops at the first zero. These are parasail's sw_trace/sg_trace choices, and pinned against them.

Output¶

Records come out in input order (@HD SO:unsorted) by default, or coordinate-sorted with --sort coordinate (see Sorted output). The header is @SQ from the reference in file order, the input's @RG, @PG and @CO copied through, and one @PG ID:eschalign for this run, chained to the last.

Every input tag is copied byte for byte — mv, ns, ts, MM/ML, RG, qs, du, ch, st, and the rest. NM, MD, AS, XS and XA are the aligner's and replace any the input carried. Of the input flags only QC-fail (0x200) and duplicate (0x400) survive.

Per read:

--secondary also writes each XA member as a 0x100 record right after its primary. Secondaries carry SEQ/QUAL as * and only NM, MD, AS and RG (minimap2's convention), so the move table and modification calls are not duplicated once per tie.

A read whose S is below --min-score — or that aligns no base at all, or is longer than --max-read-len — is written unmapped (flag 4) with its tags intact. Reads are never dropped (except by --read-ids).

--strand both also scores the reverse complement of every read. A reverse winner gets flag 16 with SEQ reverse-complemented and QUAL reversed, as SAM requires; per-base tags (mv, MM/ML) are copied verbatim, as bwa -C does. The tie order is reference order, forward before reverse.

Sorted output¶

--sort coordinate writes @HD SO:coordinate and the records in exactly samtools sort's order: by reference (@SQ order), then position, then forward before reverse strand, then input order — so the output is record for record what samtools sort of the unsorted output gives (tested on the fixtures, and checked on a 1.1 M-read sample). Unmapped reads come last, in input order; --secondary records sort by their own coordinates. Nothing else in the header changes but the @PG CL. samtools index takes the file as it is; no .bai is written.

It is a bucket sort — one bucket per reference, which is cheap because the panel is small — so there is no separate sort pass and no merge. A sorted file cannot be written until the last read is aligned, so the records are held until then; with -o - the output starts once the input ends. Past --sort-memory (default 4 GiB) the largest buckets are appended to one temporary file in --tmp-dir, and at the end each reference's records are read back and ordered on their own: peak memory is about the budget plus the largest single reference's records, not the sample. The temporary file is unlinked as soon as it is created, so nothing is left behind on any exit.

MD and NM at ambiguity codes¶

The scorer treats a reference N as matching anything. MD and NM do not: they are exactly what samtools calmd writes for the same CIGAR, because the windowed charging classifier in escapepod-rs rebuilt each read's reference from MD, and a reference N that vanished from MD would silently put the read's own base in its place. Measured on samtools 1.23.1 and pinned by a golden:

• A base is a match only when both letters are A/C/G/T (U counting as T, case ignored) and equal. A reference N — or any other IUPAC code — is a mismatch against every read base, N included, and counts toward NM.
• A mismatched or deleted reference letter is written upper-cased but otherwise verbatim (n → N; a reference U stays U).
• NM = mismatches + inserted bases + deleted bases.

So a read that matches its reference everywhere except across the reference's 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.

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:
4. both mates align at least one base;
5. 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;
6. 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.
7. 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.

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 on AVX2, 32 on AVX-512BW, i16 arithmetic that is exact (not approximate) for any panel in the intended range — chosen at runtime from what the CPU supports. The winners' tracebacks run on a second pair of kernels whose lanes each carry their own read and reference, so a read that ties with the whole panel (an adapter-only read on an adapter-flanked panel does) costs about what scoring it did. Every kernel is tested for equality against a scalar reference implementation, which is itself tested against parasail.

ESCHALIGN_BACKEND=scalar|avx2|avx512 caps the kernel choice, for A/B measurements; the run's second log line names the kernel in use. ESCAPEPOD_ALIGN_BACKEND is the same switch under its old escapepod-rs name, kept as a deprecated alias; ESCHALIGN_BACKEND wins when both are set.

Long reads¶

A tRNA sample's median read is ~130 nt, but a real one carries a handful of reads of 50–400 kb, and each costs as much as thousands of ordinary reads.

By default they are not aligned at all. --max-read-len (default 1000) sets the longest read that is scored: nothing longer is a tRNA read in any practical sense. A longer read is neither scored (on the CPU or the GPU) nor traced, but it is still written, in input order, as an unmapped record (flag 4) with every input tag copied and no NM/MD/AS/XS/XA — exactly as a read below --min-score is. The run's last lines say how many were skipped. Under --strand both neither strand is scored.

--max-read-len 0 removes the limit and aligns every read. Then two further things keep the long ones from holding up the run:

• Reads of 1,024 nt or more are scored by a second layout of the same SIMD kernel, which walks the read row by row across the references instead of reference column by column down the read. Its working set is one row of the panel, so it stays in cache at any read length, where the default layout's grows with the read — 2.3× faster on a 395 kb read, and never slower from ~1 kb up. Below 1,024 nt nothing changes.
• A read of 16 kb or more (which already gets a unit of work of its own) has its panel scored across all threads, one group of references per thread, instead of on one — so the output writer, which must emit records in input order, is not left waiting on it. With --max-read-len 0, input with many reads of 16 kb or more close together can still stall the writer.

Both give the same scores as before, so the output is identical. With --device gpu, reads over 4,096 nt are left to the CPU and take the same path. ESCHALIGN_ROW_MAJOR_ONLY=1 turns both off — every read on the original layout and one thread per unit of work — for A/B measurements only. (ESCAPEPOD_ALIGN_ROW_MAJOR_ONLY is a deprecated alias; the new name wins whenever it is set, even to 0.)

Performance¶

A 3.3 M-read dorado uBAM against the 164-reference sacCer3 dual-adapter panel aligns in ~2 minutes on 16 cores (-t 32, ~28 k reads/s, ~3,300 CPU-seconds) in under 3 GiB — against a 12 h / 160 GB budget for the bwa step it replaces. On Cascade Lake the AVX-512 and AVX2 kernels are within a few percent of each other; scalar is ~32× the CPU. Details and reproduction in benchmarks/README.md.

GPU scoring¶

In a binary built with --features gpu, scoring every read against every reference can run on a CUDA GPU: one GPU thread per (read, reference) pair, a whole --batch-size batch per kernel launch. Only the scores come from the device — the tie set, the tracebacks, MD/NM and the records are the CPU path's own code — so --device gpu and --device cpu write the same records byte for byte (tested, and checked on a 490 k-read sample). Reads over 4,096 nt are scored on the CPU inside the same run, and a panel with a reference over 3,072 nt cannot use the GPU (auto then stays on the CPU and says so).

On a gpu node (-c 16, one A30), the 3.3 M-read sample takes ~70 s and ~700 CPU-seconds with GPU scoring against 184 s and ~2,900 CPU-seconds on the same 16 cores without it, so --device auto uses the GPU when one is visible. Nothing CUDA is needed to build; at run time the GPU path needs the CUDA driver and NVRTC (pixi run install-gpu sets up the NVRTC half). Details in benchmarks/README.md.

Replacing bwa in a pipeline¶

The aa-tRNA-seq alignment step today is:

samtools view -u -N ids.txt reads.ubam \
  | samtools fastq -T '*' - \
  | bwa mem -C -H rg.sam -x ont2d -k 6 -W 13 -T 20 ref.fa - \
  | samtools view -u -F 2324 - \
  | samtools sort -o aln.bam
samtools calmd -Qb aln.bam ref.fa > aln.calmd.bam

eschalign does all of it in one process — the tags ride through without the FASTQ-comment round trip, the uBAM's @RG lines are copied without the -H header hack, MD is written at alignment time, and --sort coordinate replaces the sort:

eschalign reads.ubam -r ref.fa -o aln.bam --read-ids ids.txt --scoring 1,-1,-2,-1 \
    --sort coordinate
samtools index aln.bam

escpod classify (escapepod-rs) reads either order. Unmapped reads are written, never dropped; samtools view -F 2324 reproduces the old output's contents exactly (primary, forward, mapped records) if that is wanted, and a filtered sorted file stays sorted. Stamping sample/barcode read groups stays pipeline plumbing.

The scoring scheme is the pipeline's call, judged by charging-call agreement: --scoring 1,-1,-2,-1 is bwa -x ont2d's penalties, but bwa also clips, bands and chains, so the same scores do not give the same CIGAR for every read. On escapepod-rs's classify test fixture (the same 60 reads as crates/eschalign-cli/tests/fixtures), 27 of 60 reads get a different CIGAR (or start) from bwa's; the charging calls of every read whose alignment is unchanged are bit-identical.

Examples¶

# dorado uBAM in, aligned BAM out
eschalign calls.ubam -r trna.fa -o aligned.bam -t 32

# FASTQ, overlap mode, gap-lenient scores
eschalign reads.fq.gz -r trna.fa -o aligned.bam --mode semiglobal --scoring 1,-1,-2,-1

# Coordinate-sorted, ready for samtools index
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¶

Deprecated aliases from escapepod-rs (escpod align): ESCAPEPOD_ALIGN_BACKEND and ESCAPEPOD_ALIGN_ROW_MAJOR_ONLY. They still work; the ESCHALIGN_* name wins whenever both are set.

Logging¶

Status goes to stderr through tracing — info by default, -v debug, -vv trace, -q errors only; dependencies are held at warn. Only the BAM (under -o -) goes to stdout.