Skip to content

[auto-optimize]: Speed up alignment about 3x with byte-identical output - #536

Open
jmschrei wants to merge 3 commits into
BenLangmead:masterfrom
jmschrei:speed-up-alignment
Open

jmschrei wants to merge 3 commits into
BenLangmead:masterfrom
jmschrei:speed-up-alignment

Conversation

@jmschrei

Copy link
Copy Markdown

This optimization was found with auto-optimize, a tool developed by the Programmable Genomics Laboratory at UMass Chan Medical School: it speeds up bowtie2-align by 1.7x to 4.3x on every input timed except a 12-read one at -p 1 (1.0x, where index loading dominates), with output byte-identical to master.

On the four main cases at -p 4 (ATAC-seq, ChIP-seq, CUT&RUN and DNase-seq against hg38), the mean wall time falls from 19.70 s to 6.49 s (3.04x). Most of the gain comes from three things: rewriting the 8-bit dynamic-programming fills, batching and overlapping FM-index walks across reads, and loading the index in parallel. Peak memory is unchanged apart from small --mm runs (see Memory), and no OS-specific calls are added.

It builds on #535, which adds a pytest suite and changes no source. The first commit here is #535's. This pull request's own change is the second commit: 37 files, +7,472 / -1,815. Once #535 is merged, only the second commit remains.

Changes

In the order they were made:

  1. End-to-end 8-bit fill (aligner_swsse_ee_u8.cpp): an exact prefix scan across lanes replaces the lazy-F loop, which took about 18 passes per column.
  2. AVX2 dispatch (bowtie_main.cpp, Makefile): the dispatch tested AVX2 in CPUID leaf 1 instead of leaf 7, so bowtie2-align-*-v256 ran only on CPUs where an unrelated bit happened to be set. It now requires every x86-64-v3 feature (CPUID leaves 1, 7 and 0x80000001) and that the OS saves the AVX registers (XCR0). The -v256 build compiles with -ffp-contract=off, so its floating-point results, and with them --score-min and MAPQ, match the 128-bit build. Short reads get a 128-bit fill path, where the 256-bit fill was slower.
  3. Local-mode 8-bit fill: the candidate cells are gathered with vector instructions after the fill.
  4. One-mismatch seed searches: 16 independent FM-index walks are stepped in turn, each with a prefetch of its next BWT side, so their memory waits overlap.
  5. End-to-end 8-bit fill: a shorter dependency chain from one column to the next, 107 to 41 cycles per column in the inner loop.
  6. Local-mode 8-bit fill: the same chain shortening, and the lazy-F loop removed.
  7. Local-mode candidate processing: a counting sort instead of a comparison sort, and a vectorized test of whether a candidate is dominated.
  8. Seed-hit cache: when it records a BW range, its list of unresolved offsets is filled a page-sized run at a time (std::fill_n, new PList::addFillUpTo in ds.h) instead of with one add() call per element.
  9. Index arrays and short fills:
    • The large index arrays are allocated 64-byte aligned, so a BWT side read at a random offset touches one cache line instead of two.
    • The end-to-end fill for reads of up to 32 rows runs as two interleaved dependency chains.
  10. Mate sweeps: the exact end-to-end sweeps of both mates of a pair are interleaved.
  11. Query profile: the striped query profile is built by table lookup.
  12. 8-bit backtraces: they clear only the intervals they wrote, instead of zeroing whole matrices. master zeroes 14.4 GB per 40,000 ATAC-seq pairs.
  13. Exact sweeps across reads: the exact end-to-end sweeps are batched across the reads of an input batch, so walks of different reads overlap.
  14. SA resolution, single rows: single-row walks run ahead and are committed where the serial walk would use them.
  15. SA resolution, large ranges: walks of large BW ranges proceed one 64-bit BWT word at a time.
  16. SA resolution, steps inside range walks: the single-row steps inside range walks also run ahead.
  17. Occurrence counts: branch-free in every LF step.
  18. RedundantAlns: keeps each read's cells as runs along diagonals instead of a set of cells.
  19. Concurrent index loading (new index_load.h): at -p > 1, the four large index arrays load concurrently on up to four helper threads. File errors are reported at the same point and with the same messages as the serial load.
  20. BtBranchTracer: resets lazily, only when the checkpointed backtracer runs.
  21. One-mismatch searches across reads: batched across the reads of an input batch.
  22. Reference windows: unpacked from the 2-bit reference directly into the masks the aligner uses.
  23. Prefetch timing: ftab lookups are prefetched one round early, and a sweep restart's first LF step moves to the next turn.
  24. Seed-hit ranking and sampling: without linear scans or a sort per call. The random draws are unchanged.
  25. Two-chain fill, longer reads: the two-chain end-to-end fill extends to reads of 33-128 rows, which covers mate rescue in ATAC-seq.
  26. Sampled SA elements: built on first use (94-97% are never used). The random draws are unchanged.
  27. Both strands together: the seed searches of both strands of a read run in one batch, and per-call scratch space is kept between calls.
  28. Seed instantiation: redone only for inputs that changed (98.5% of seeds skip it).
  29. Seed-hit extension: its left and right LF chains are stepped in turn.
  30. joinedToTextOff: uses a bucket table built at index load instead of a binary search over the reference fragments.
  31. Local-mode backtrace: runs without a frame stack, and the domination test reads a sentinel-padded array.
  32. Output: the output queue no longer writes shared state for every read, SAM fields of known length are written directly, and reverse and reverse-complement are fused.
  33. Batched one-mismatch prefetches: each round's BWT-side prefetches are issued together, after the round's steps.
  34. Index loading in pieces: the arrays load in 64 MB pieces from one queue, and the main thread reads pieces too. The error messages of a failed read are unchanged.

Tests and CI

The pytest suite from #535 passes on this branch: 611 passed, 1 skipped. The Linux jobs of the repository's two workflows were run locally on this commit, as the workflow files write them, each in a fresh clone:

Job Result
make clean; make -j4 allall; make simple-test passed
the same with SSE_AVX2=0 passed
make clean; make -j4 allall; NUM_CORES=1 make random-test passed
the same with SSE_AVX2=0 passed
python3 -m pytest scripts/test on both builds 611 passed, 1 skipped

The macOS jobs could not be run here.

Behaviour changes

There is one. On CPUs where the old dispatch failed to detect AVX2, bowtie2-align-s and -l now run their -v256 builds, and with -ffp-contract=off those give the same output as the 128-bit builds. On CPUs where the old dispatch did select -v256, fused multiply-adds could change MAPQ relative to the 128-bit build. We measured that in 7,886 of 200,000 24-bp reads with a custom --score-min. Users on such CPUs therefore now get the 128-bit build's output. With default settings we saw no such difference: 221 of 221 jobs were byte-identical against master's own -v256 binary.

Usage benchmarks

Machine: AMD EPYC 9575F, 2 sockets x 64 cores (256 hardware threads, 2 NUMA nodes), 1 TB RAM, Ubuntu 24.04.5, Linux 7.0, GCC 13.3.0, Perl 5.38.2. Both builds use the repository's make.

Inputs: no synthetic data. Every input is public sequencing data, aligned against hg38 (UCSC, 455 sequences) indexed by master's bowtie2-build.

  • Main cases:
    • ATAC-seq HCT116 (ENCODE ENCSR260SWI).
    • ChIP-seq K562 CTCF (ENCSR000EGM), converted from Phred+64.
    • CUT&RUN K562 CTCF (GEO GSE126612, SRR8581590), with the Henikoff settings --end-to-end --very-sensitive --no-mixed --no-discordant -I 10 -X 700.
    • DNase-seq T cells (ENCSR025EYJ).
  • Other assays: 250,000-read samples from public ENCODE, 4DN, GEO/SRA and 10x Genomics data.

How it was timed:

  • Isolation: each process ran alone, pinned to one NUMA node (numactl --cpunodebind=0 --membind=0), and started only when at most two processes were runnable on the whole machine.
  • Repeats: each configuration ran 3 times per build (once for the two full-depth inputs), alternating which build went first. The table gives the median wall time.
  • What the time includes: index and reads were in the page cache, and SAM went to local disk. Times cover wrapper start-up, index loading, alignment and output.
  • Peak memory: the process's maximum resident set (/usr/bin/time %M).
  • Output check: whether the SAM was byte-identical, after removing @PG and sorting records when -p > 1 without --reorder.
Input class Input Threads master (s) This branch (s) Speedup master peak (MB) This branch peak (MB) Output
main cases ATAC-seq HCT116, PE101, 400k pairs, -X 2000 1 100.18 25.71 3.90x 3,418 3,418 identical
main cases ATAC-seq HCT116, PE101, 400k pairs, -X 2000 2 51.21 14.05 3.64x 3,436 3,434 identical
main cases ATAC-seq HCT116, PE101, 400k pairs, -X 2000 4 26.98 7.91 3.41x 3,466 3,470 identical
main cases ATAC-seq HCT116, PE101, 400k pairs, -X 2000 8 14.79 4.84 3.06x 3,534 3,532 identical
main cases ATAC-seq HCT116, PE101, 400k pairs, -X 2000 16 8.61 3.37 2.55x 3,660 3,658 identical
main cases ChIP-seq K562 CTCF, SE36, 678k reads 1 24.43 10.12 2.41x 3,410 3,412 identical
main cases ChIP-seq K562 CTCF, SE36, 678k reads 2 13.03 5.58 2.34x 3,412 3,416 identical
main cases ChIP-seq K562 CTCF, SE36, 678k reads 4 7.45 3.29 2.26x 3,428 3,434 identical
main cases ChIP-seq K562 CTCF, SE36, 678k reads 8 4.52 2.20 2.05x 3,454 3,458 identical
main cases ChIP-seq K562 CTCF, SE36, 678k reads 16 3.05 1.58 1.93x 3,500 3,506 identical
main cases CUT&RUN K562 CTCF, PE25, 250k pairs, Henikoff settings 1 81.93 19.02 4.31x 3,416 3,418 identical
main cases CUT&RUN K562 CTCF, PE25, 250k pairs, Henikoff settings 2 41.70 9.97 4.18x 3,430 3,432 identical
main cases CUT&RUN K562 CTCF, PE25, 250k pairs, Henikoff settings 4 21.81 5.39 4.05x 3,460 3,456 identical
main cases CUT&RUN K562 CTCF, PE25, 250k pairs, Henikoff settings 8 11.58 3.14 3.69x 3,518 3,512 identical
main cases CUT&RUN K562 CTCF, PE25, 250k pairs, Henikoff settings 16 6.52 2.07 3.15x 3,622 3,614 identical
main cases DNase-seq T cell, PE76, 401k pairs, --local 1 81.22 32.18 2.52x 3,420 3,420 identical
main cases DNase-seq T cell, PE76, 401k pairs, --local 2 41.83 17.07 2.45x 3,438 3,440 identical
main cases DNase-seq T cell, PE76, 401k pairs, --local 4 22.57 9.36 2.41x 3,474 3,474 identical
main cases DNase-seq T cell, PE76, 401k pairs, --local 8 11.72 5.59 2.10x 3,540 3,544 identical
main cases DNase-seq T cell, PE76, 401k pairs, --local 16 7.29 3.64 2.00x 3,666 3,694 identical
other assays ATAC-seq A549, SE51, 250,000 reads 4 4.20 1.77 2.37x 3,434 3,434 identical
other assays CAGE GM12878, SE36, 250,000 reads 4 2.88 1.29 2.23x 3,424 3,428 identical
other assays ChIP-seq BLaER1 H3K4me1, SE50, 250,000 reads 4 4.36 1.76 2.48x 3,434 3,438 identical
other assays CUT&RUN HUVEC BRG1, PE36, 250,000 pairs 4 12.85 4.19 3.07x 3,452 3,450 identical
other assays CUT&Tag K562 H3K27me3, PE25, 250,000 pairs 4 16.38 4.40 3.72x 3,442 3,442 identical
other assays DNase-seq GM12878, SE36, 250,000 reads 4 3.17 1.47 2.16x 3,428 3,430 identical
other assays eCLIP K562 RBFOX2, PE45/46, 250,000 pairs 4 11.74 3.93 2.99x 3,452 3,452 identical
other assays Hi-C K562, PE101, 250,000 pairs 4 23.93 6.92 3.46x 3,468 3,462 identical
other assays MNase-seq K562, PE50, 250,000 pairs 4 14.47 4.52 3.20x 3,458 3,454 identical
other assays 10x Multiome ATAC, PBMC, 250,000 pairs 4 10.89 3.67 2.97x 3,456 3,452 identical
other assays PRO-cap HCT116, PE150, 250,000 pairs 4 16.93 6.64 2.55x 3,460 3,464 identical
other assays PRO-seq K562, SE, 250,000 reads 4 7.88 2.72 2.90x 3,436 3,438 identical
other assays Ribo-seq K562, SE51, 250,000 reads 4 3.31 1.29 2.57x 3,428 3,430 identical
other assays RNA-seq K562 poly(A), SE76, 250,000 reads 4 4.26 1.92 2.22x 3,432 3,436 identical
other assays STARR-seq K562 RNA, PE100, 250,000 pairs 4 22.00 6.58 3.34x 3,468 3,468 identical
other assays WGS HG002, PE148, 250,023 pairs 4 31.79 12.50 2.54x 3,482 3,474 identical
presets and options ATAC-seq HCT116, --very-fast 4 22.63 6.76 3.35x 3,462 3,462 identical
presets and options ATAC-seq HCT116, --very-sensitive 4 37.98 10.78 3.52x 3,496 3,494 identical
presets and options ATAC-seq HCT116, --local 4 31.45 13.26 2.37x 3,488 3,484 identical
presets and options ATAC-seq HCT116, --very-sensitive-local 4 40.13 16.53 2.43x 3,498 3,502 identical
presets and options ATAC-seq HCT116, -k 5 4 37.67 9.99 3.77x 3,472 3,470 identical
presets and options ATAC-seq HCT116, --reorder 4 27.45 8.03 3.42x 3,468 3,470 identical
presets and options ChIP-seq K562 CTCF, -k 10 4 9.16 5.16 1.78x 3,434 3,432 identical
presets and options ChIP-seq K562 CTCF, --local 4 8.23 3.73 2.21x 3,442 3,440 identical
presets and options ChIP-seq K562 CTCF, --mm 4 8.72 4.12 2.12x 3,358 3,356 identical
presets and options ChIP-seq K562 CTCF, --very-sensitive 4 8.88 3.95 2.25x 3,442 3,444 identical
tiny input ChIP-seq BLaER1 H3K4me1, 12 reads 1 0.86 0.84 1.02x 3,406 3,404 identical
tiny input ChIP-seq BLaER1 H3K4me1, 12 reads 4 0.91 0.47 1.94x 3,406 3,406 identical
tiny input ChIP-seq BLaER1 H3K4me1, 12 reads 16 1.05 0.61 1.72x 3,408 3,410 identical
full depth ATAC-seq HCT116, PE101, 47,166,588 pairs 16 859.15 218.75 3.93x 3,748 3,766 identical
full depth ChIP-seq K562 CTCF, SE36, Phred+64, 33,894,235 reads 16 79.67 33.96 2.35x 3,565 3,518 identical
The SAM was byte-identical in every configuration.

On the 105 full-depth experiments below, the summed wall time was 25.4 h for master and 7.9 h for this branch (3.20x). Those runs went two or three at a time beside other jobs, so they are not a controlled timing.

How the output was checked

The usage benchmarks above check speed and memory. Separately, we ran master and this branch side by side on a collection of public sequencing data and hand-made malformed files:

  • Scale: 7,109 distinct command lines, run under up to four pairs of builds. That is 18,087 comparisons: 36,876 process runs and 123 hours of process time, including the controls.
  • What was compared: for every command, the exit code, stdout, stderr and every file written, byte for byte. Build and output paths were replaced by tokens first, so @PG lines are compared too.
  • Multithreaded order: when -p > 1 without --reorder, record order varies from run to run in master itself, so those runs were compared as sorted records.

Builds compared

Pair master This branch Why
Release make make what users run. On this CPU, master runs its 128-bit code (the dispatch bug) and this branch its -v256 code
SSE2-only make make SSE_AVX2=0 this branch's 128-bit code
ASan + UBSan bowtie2-align-*-sanitized, built with -fsanitize-recover same memory errors and undefined behaviour
TSan -fsanitize=thread same data races
master -v256 bowtie2-align-*-v256 run directly make master's AVX2 code, which its dispatch never selects here
Controls make master again whether a difference also appears between two runs of master

What ran

Suite Commands What
Read tiers 833 every human FASTQ in the collection at tiny to large size (5 to 1,000,601 reads per file; 21 assays, single- and paired-end), default settings, -p 1
--local, --very-sensitive, --phred64 559 the same files at one size each, under these settings; --phred64 on the Phred+64 files
Threads 108 -p 3, 4, 8 and 16, with and without --reorder
Options 296 input formats (FASTQ plain, gz and bz2; FASTA; raw; -c; -F; tab5/tab6; qseq; interleaved; stdin; comma-separated lists; mixed paired and unpaired; BAM, aligned and unaligned, with --align-paired-reads and --preserve-tags), trimming and read selection, presets, scoring, reporting (-k, -a), pairing, outputs (--un*/--al*, --met-file, SAM options), threads, indexes, and command-line errors
Flag values 1,760 each of the 262 options of bowtie2, bowtie2-build and bowtie2-inspect with several values each, single-end, paired-end and -p 4 --reorder
Combinations 514 a pairwise covering array over option values: every pair of values of two options occurs in some command
bowtie2-build, bowtie2-inspect 640 builds of 12 genomes (phage lambda to mm10) and of every malformed FASTA file, each builder option, and inspect on the results
Malformed reads 2,036 every malformed and edge-case FASTQ (513), BAM/SAM (628) and FASTA (895) in the collection, e.g. truncated or empty files, CRLF line ends, missing + lines, unequal mate counts, zero-length reads, lowercase bases, odd quality characters, gzip CRC errors, BGZF without its EOF block
Damaged indexes 246 80 damaged copies of the lambda and hg38 indexes (files truncated, emptied, missing, padded or bit-flipped; forward and mirror files swapped; a large index renamed; reference files from another genome) at -p 1, -p 4 and -p 4 --mm
Long reads 12 ONT and PacBio reads, under a time limit
Full depth 105 every full-depth experiment in the collection except 10x barcode reads: 9.76 billion reads, 20 assays, hg38 (one mouse and one yeast experiment against mm10 and sacCer3), -p 16 --reorder
TSan 70 -p 4 and -p 16, ordered and unordered output, damaged indexes
master -v256 221 the small read tier (220 commands) and one full-depth DNase-seq experiment (27M pairs, --local)

Results

Pair Commands Identical Identical failure Differ
Release 7,109 5,833 993 282
SSE2-only 5,385 4,305 1,001 79
ASan + UBSan 5,523 3,174 455 1,894
TSan 70 0 6 64
master -v256 221 221 0 0
Controls (master vs master) 129 7 9 113

"Identical failure" means both builds failed with byte-identical output. None of the differences comes from a change in what bowtie2 computes. Each cause is listed below with the evidence for it:

Cause Evidence
Output that changes from run to run by design: warnings from threads at -p > 1, multithreaded bowtie2-build logs, --met-file, --met-stderr, --read-times, --non-deterministic, multithreaded --log-dp master differs from itself the same way (113 of the 129 controls)
@PG command lines that record the wrapper's temporary file name (it carries the wrapper's PID, e.g. for bz2 input) or the output path of a --bam run the difference is confined to @PG
--version, which prints the compiler options of the binary that runs 3 commands
bowtie2-inspect on an index built with -r prints an uninitialized value 20 commands; both builds read the uninitialized value
An assertion message with a different source line number (bt2_idx.h:211 vs :297) 2 commands
Both builds crash, abort or time out on malformed input or a bug in master; the records written before that are identical 153 release commands; same exit code in every one
Under ASan, both builds hit a 24 GB per-job memory cap on 31 kb FASTA records in --local mode 368 commands
Sanitizer reports (below) 1,419 ASan/UBSan and 62 TSan commands that are identical once the reports are removed

Apart from these, every exit code and output file is the same on valid input. On malformed input, both builds fail the same way in every case but one. When an index's .3/.4 files come from a different genome than its .1/.2 files, master segfaults and the SSE2-only build of this branch loops. Behaviour there is undefined in both builds, because nothing checks that the files belong together.

Sanitizers

  • Only in this branch: none.
  • Only in master, and no longer triggered by this branch:
    • a stack-buffer-overflow and an out-of-bounds index in SeedAligner::oneMmSearch on reads with N (aligner_seed.cpp, 545 commands);
    • misaligned and short loads in Ebwt::countUpTo (556 commands).
  • In both builds, in the same commands:
    • a null source pointer passed to memcpy in SStringExpandable::install;
    • misaligned TIndexOffU loads with --mm on a large (.bt2l) index;
    • a heap overflow on a FASTQ record without its + line;
    • overflows reading --preserve-tags BAM input;
    • leaks, shifts and divisions by zero on damaged indexes.
  • TSan: one data race, the same in both builds. OutFileBuf::flushBlocking() reads asyncData_.buf without the lock that AsyncData::writeComplete() holds when it clears it.

Threads at -p 1

The speedup does not come from extra threads:

  • Thread count: at -p 1, both builds create the same three threads, one worker plus the input reader and output writer that master already runs. strace counted them on single-end ChIP-seq and on paired-end --local DNase-seq.
  • CPU time: both builds use 1.00 to 1.04 CPU-seconds per wall second, and this branch's CPU time falls with its wall time. On ATAC-seq at -p 1, for example, it drops from 101.5 to 26.8 CPU-seconds.
  • At -p > 1: this branch adds up to four short-lived threads that read the index files in parallel. With -p 1, --mm, --shmem, --verbose, -t or --sanity, the index loads serially as before.

Memory

Peak resident memory is unchanged: 3.3 to 3.8 GB for hg38 in every timed configuration, a ratio of this branch to master between 0.987 and 1.007.

Across the 4,390 identical release commands that peaked above 500 MB without --mm, the median ratio is 1.000 and the largest is 1.018. Seven used less than 95% of master's memory, and all seven are -a runs.

With --mm, where the index is memory-mapped instead of copied, the resident set also counts the index pages a run touches. Over the 102 --mm commands the ratio ranged from 0.41 to 2.17:

  • The most extreme case: 1,509 CUT&RUN pairs with --local --mm -p 3. This branch peaked at 1.20 GB against 0.55 GB in three repeats, and took 0.31 s against 0.19 s.
  • What the extra memory is: file-backed and shared (RssFile 947-1,044 MB against 265-346 MB, sampled every 20 ms). Private memory was the same, 17-18 MB.
  • A larger input: on the 677,596-read ChIP-seq input with --mm, both builds reach 3.3 GB, and this branch is 2.12x faster (table above).

We did not establish which change touches more of the mapped index at start-up.

What this testing does not cover

  • Platforms: only x86-64 Linux, one AMD Zen 5 CPU and GCC 13. We tested no Intel CPU, ARM, macOS, Windows or other compiler.
  • Not compared: two small long-read inputs, which reach the 30-minute limit in every build; one bowtie2-build option combination, which was too slow; four ASan jobs that hang in both builds (-F on zstd and xz FASTA).
  • Non-default --score-min on master's -v256 code: see Behaviour changes. Only the default-settings check above was run.

Bugs in master found along the way

We reproduced 25 bugs in master, and this branch behaves the same on each. We will file them as separate issues. Among them:

  • -F k:<int>,i:<int>, the syntax --help documents, yields 0 reads.
  • --int-quals aborts on any multi-base read.
  • bowtie2-build --nodc and --threads 0 die with SIGFPE, while the wrapper exits 0.
  • Paired-end plus unpaired input with --reorder at -p > 1 hangs.
  • With -b, the --un/--al files contain raw BAM records.
  • The stale MC:Z tag makes output nondeterministic at -p > 1.

🤖 Generated with Claude Code

Jacob Schreiber and others added 2 commits October 16, 2026 10:32
… CLIs

612 tests in scripts/test/test_*.py covering input formats and read
preprocessing, alignment, scoring, reporting and effort options,
paired-end semantics, output and SAM options, the wrapper scripts, and
index building and inspection. Expected values are derived from reads
cut at known positions from fixed random references, using the MANUAL's
scoring rules; tests that pin current output say so in their docstring.

Run after building, from the repository root:

	make -j4 allall
	python3 -m pytest scripts/test

Optional tools (zstd, lz4, bzip2, samtools) and the -debug/-sanitized
binaries are skipped when absent. BOWTIE2_BIN_DIR points the suite at
binaries built elsewhere.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
bowtie2-align is 1.7x to 4.3x faster on every input timed (end-to-end
and local, single- and paired-end, -p 1 to -p 16), except a 12-read
input at -p 1, where index loading dominates. Its output is byte for byte
the same as before. Peak memory is unchanged, except that with --mm a
small run touches more of the mapped index.

Dynamic programming (aligner_swsse_ee_u8.cpp, aligner_swsse_loc_u8.cpp,
aligner_swsse.cpp, aligner_sw*.cpp):
- The 8-bit end-to-end fill computes the vertical-gap term with an exact
  prefix scan across lanes instead of the lazy-F loop. Both 8-bit fills
  have a shorter dependency chain from one column to the next. Reads of
  up to 128 rows are filled as two interleaved chains.
- The query profile is built by table lookup. The 8-bit backtraces clear
  only the intervals they wrote, not whole matrices.
- Local mode gathers its candidate cells with vector instructions, sorts
  them by counting, tests domination vectorized, and backtraces without
  a frame stack.

FM index and seeds (aligner_seed.cpp, group_walk.h, bt2_idx.*):
- Occurrence counts in every LF step are branch-free, and the large index
  arrays are 64-byte aligned.
- Exact and one-mismatch seed searches are batched across the reads of an
  input batch and across both strands. Independent walks are stepped in
  turn, with each walk's next BWT side prefetched, so memory waits overlap.
- Seed instantiation is redone only for inputs that changed. Seed-hit
  ranking and sampling no longer scan lists linearly or sort per call.
- SA resolution runs single-row walks ahead of their use, and walks large
  BW ranges a 64-bit BWT word at a time. Sampled SA elements are built on
  first use. Every random draw is unchanged.
- joinedToTextOff uses a bucket table built at index load instead of a
  binary search over the reference fragments.

Elsewhere: RedundantAlns keeps per-read runs along diagonals,
BtBranchTracer resets lazily, reference windows are unpacked straight to
masks, and the output queue and SAM writer do less work per record.

Index loading (index_load.h): at -p > 1, the four large arrays (.4 file,
forward and mirror BWT, SA sample) are read in 64 MB pieces by up to four
helper threads and the main thread. Read errors print the same messages as
the serial load. The load stays serial at -p 1, and with --mm, --shmem,
--verbose, -t or --sanity.

AVX2 dispatch (bowtie_main.cpp, Makefile): the dispatch tested AVX2 in
CPUID leaf 1, so bowtie2-align-*-v256 ran only on CPUs where an unrelated
bit was set. It now requires every x86-64-v3 feature (CPUID leaves 1, 7
and 0x80000001) and OS support for the AVX registers (XCR0). The -v256
build compiles with -ffp-contract=off, so its floating-point results,
--score-min and MAPQ, are rounded as in the 128-bit build.

No OS-specific calls are added.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@jmschrei

Copy link
Copy Markdown
Author

Hi Ben, hope you're doing well. The above message was written by Claude, but this comment is written by me (though based on my MACS3 PR).

I have been developing a tool, called [auto-optimize], which basically applies the auto-research loop to code to improve speed and memory usage of code, but with extensive robustness checks to try to ensure nothing subtly breaks. I have been applying [auto-optimize] across the genomics ecosystem, starting with my own tools but now extending to other important ecosystem tools, with the hope of making things broadly better. Most recently, an auto-optimization of MACS3 was merged into that repository: macs3-project/MACS#756, and I've been impressed with other findings, such as this 20 speedup on my FIMO implementation: jmschrei/memesuite-lite#22.

In an effort to be more robust, the process first augments the unit tests to make sure that nothing breaks during the search. See #535 (feel free to merge that independently of any choice you make here). Any iteration in our loop that broke any unit test was automatically discarded, and we used a set of external validations that were not directly optimized against each iteration to make sure we got identical results. At the end, we ran several thousand mapping jobs across assays and data type (single-end vs paired-end, read length, etc), covering every command-line flag, combinations, etc. Importantly, we also created a set of malformed data (both the reads and the genome) to make sure that not only did you get identical outputs on good data, but that it failed the same way on bad data so that we didn't just get rid of safeties that were necessary but slow.

If you look at the results below, you'll see that there were dozens of small improvements that, together, led to the large improvements in speed. The big dip at the start was the end-to-end 8-bit DP fill. The purple run at the end was me reverting one change that required massive pages and may have caused issues on some architectures, which was probably a barrier to adoption.

image

In an effort to be consistent with https://rewrites.bio/, I'd like this to be merged into your repository so you can continue to get credit for all the work you've done and work you've done for the field (though I don't think that's really in question...). However, as a maintainer myself, I understand being skeptical of merging a PR of this size, particularly when the submitter can't explain everything in it (a consequence of so many changes being made at a low level). If you feel uncomfortable with it, we're happy to host a new repository that would point all credits and citations back to you but be a stand-alone entity. We're happy to share any logs generated during the process with you as well, or anything else you might need.

In the 8-bit end-to-end and local backtraces, the step from vector row 0
to the previous element of the last vector row is btUpWrap =
(colstride - ROWSTRIDE) * NBYTES_PER_REG - 1, a size_t. When the matrix
has one vector row (reads of up to 16 bases with 16-byte vectors), it is
(size_t)-1, and btcur + btUpWrap wraps around the address space to
btcur - 1. The address is the intended one, but the wrap is undefined
behaviour, which clang's -fsanitize=pointer-overflow reports; the
sanitized build in the macOS CI job aborted on it. Casting to ptrdiff_t
at the addition makes the step a plain -1. With gcc 13 on x86-64, the
release aligners' machine code is identical with and without the cast.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@jmschrei

Copy link
Copy Markdown
Author

This reply was written with Claude Code.

The macOS "Simple Tests" failure was a real bug in this PR, fixed in 9bdab09. In the rewritten 8-bit backtraces (aligner_swsse_ee_u8.cpp and aligner_swsse_loc_u8.cpp), the step from vector row 0 to the row above, btUpWrap, is a size_t. When a read fits in one vector row (up to 16 bases), it equals (size_t)-1, and btcur + btUpWrap wraps around the address space to btcur - 1. That is the intended cell, so alignments were correct, but the wrap is undefined behaviour, and clang's UBSan aborted the sanitized build on it. GCC's UBSan does not report this pattern, which is why the Linux jobs passed.

The fix casts the offset to ptrdiff_t where it is added, so the step is a plain -1:

-		uint8_t * const btup = btwrap ? btcur + btUpWrap : btcur - btUp;
+		uint8_t * const btup = btwrap ? btcur + (ptrdiff_t)btUpWrap : btcur - btUp; // (size_t)-1 when nvecrow_ is 1

With GCC 13 on x86-64, the release aligner binaries' machine code is byte-identical with and without the cast, so speed and output are unchanged. With clang 19 on Linux, make simple-test failed on the same test as the macOS job before the fix and passes after it, with the sanitizers set to stop at the first report.

🤖 Generated with Claude Code

@jw1912

jw1912 commented Oct 20, 2026

Copy link
Copy Markdown

@jmschrei why make one massive PR? You could break this up into more manageable chunks, or ask Claude to do it for you. I think, for example, all of the e2e 8-bit fill optimisations (#1, #5, #9 second point, #25) can be nicely put into their own self-contained PR and they achieve a notable portion of the speedup.

This is also important for the stacks of small improvements, as I found that a couple of them (#23 and #33) were regressions on my hardware and a few more were not measurable improvements. I just happened upon this PR so please lmk if the hardware you are testing on is particularly relevant compared to my consumer CPU, or any other misconceptions I may have.

@jmschrei

Copy link
Copy Markdown
Author

Hi @jw1912, thanks for commenting.

The output, as provided, is what was discovered by the auto-optimization algorithm through its search process. I'd be happy to split this into smaller chunks if that's what the core developers here are interested in. I'm also happy for them to cherry-pick any individual or set of ideas into cleaner PRs.

With respect to the potential for regressions: the auto-optimization search process is supposed to only consider ideas that work on all major operating systems and not to rely on OS-specific tools or flags, but admittedly that is prompt-based so has the potential for error. It is possible that some would not transfer across machines -- it is also possible that some only become improvements once other portions of the code have been improved.

In general, I'm happy to redo this process with whichever sorts of constraints would make the most sense (e.g., increasing the noise band so only big improvements get in, avoiding certain portions of the code/flags, etc). I was more interested in starting the conversation about the best way to integrate AI-discovered speedups into these code-bases, rather than saying that this is the only way to move forward.

@BenLangmead

Copy link
Copy Markdown
Owner

@ch4rr0 and I will be taking a look at this soon. When I first read through the PR, I definitely left with the impression that it should be split up, but we'll comment more thoroughly when we've had a chance to digest and discuss.

@jmschrei

Copy link
Copy Markdown
Author

Sounds good, thanks for taking the time to discuss. Let me know if I can be helpful in any manner.

@BenLangmead

Copy link
Copy Markdown
Owner

We are indeed inclined to break this up. As a sample, here's one piece: #542.

@jmschrei -- you are co-author on all but one of those commits in PR542. That last commit is a small fix to avoid an assert in the debug binary.

Let us know if you have any comments on this. I intend to try to keep going this way; it's not clear to me yet whether we'll end up using everything from #536; but we'll keep plugging away and taking chunks.

@jmschrei

jmschrei commented Oct 24, 2026 •

Copy link
Copy Markdown
Author

Sounds great! Thanks for managing the process. Take/discard whatever you'd like.

@BenLangmead

Copy link
Copy Markdown
Owner

Here's the next one: #544

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants