A History of Pairwise Alignment

Ragnar {Groot Koerkamp}

CiE 2026, Trier

My blog: curiouscoding.nl

Pairwise alignment:
finding differences between strings

Covid – \(\alpha\), December 2020

Covid – \(\omicron\), December 2021

Pairwise alignment

• Find the mutations between two sequences

Dynamic programming

Dynamic programming

Dynamic programming

Dynamic programming

Dynamic programming

50 000 TB of sequenced DNA → it must be fast

Data Center
PhonlamaiPhoto | istockphoto.com

Growth of SRA

Goal:
Fast code

Goal:
High throughput code

Goal:
Optimal throughput code

What is high troughput code?

1. Complexity

   • Few operations:

     \(\quad O(n^2)\quad\longleftrightarrow\quad O(n)\)

2. Efficiency

   • Fast operations:

     memory read, 100 ns \(\quad\longleftrightarrow\quad\) 0.1 ns, addition

3. Implementation

   • Parallel operations:

     SIMD, instruction-level parallelism

Lots of work on this problem!

Needleman-Wunsch – Quadratic \(O(n^2)\)

Can we do better than quadratic?

No! (not much anyway)

• There is a \(O(n^{2-\varepsilon})\) lower bound assuming SETH.
• One of the best algorithms is four Russians in \(O(n^2 / \lg n)\):

Can we do better than quadratic? – take 2

Yes! BFS / Dijkstra – \(O(ns)\) for distance \(s\)

• Visit only states at distance \(\leq s\), within the central \(2s+1\) diagonals.

Diagonal transition – \(O(n + s^2)\)

• Greedily extend matches along diagonals → visit \({\leq} s\) states per diagonal

Diagonal transition – \(O(n + s^2)\)

A*PA – near-linear, \(3\times\) faster on similar seqs

[1] Exact Global Alignment Using A* with Chaining Seed Heuristic and Match Pruning.
RGK and Pesho Ivanov, Bioinformatics 2024.

A*PA heuristics

A*PA – not quite linear, and terrible efficiency

Intermezzo: array indexing is not \(O(1)\)!

Band Doubling – \(O(ns)\)

• Compute all states within central \(t\) diagonals for exponentially growing \(t\). Super efficient!

Optimizations

Bitpacking (Myers' 99)

• Computes \(w=64\) states at once by encoding the \(\{-1, 0, +1\}\) differences between adjacent states in a word.
• \(O(ns/w)\)

SIMD: single instruction multiple data, special CPU instructions.

• Operates on e.g. 4 or 8 64-bit words at the same time..

A*PA2 – good efficiency: up to \(19\times\) faster

• Band-doubling + bitpacking + SIMD + A*PA

[2] /A*PA2: Up to 19x Faster Exact Global Alignment. WABI 2024.

A*PA: comparison

But: Long-read global alignment is just not a problem people have…

Alignment modes

Semi-global variants

… Sassy: 100bp semi-global alignment

Approximate String Matching

• Find pattern \(P\) in text \(T\).
  • \(m := |P|\), \(n := |T|\)
• Allow up to \(k\) unit-cost errors.
• Find all matches.

• IUPAC support: Handle ACTG, N, YR...

• Not semi-global alignment
• Not \(\Theta(nm/w)\)
• Not new: Lots and LOTS of old literature

ASM in \(O(\lceil n/w\rceil k)\)

• Myers' bitpacking on \(w=64\)-bit blocks.
• "Early break" when all states in block have cost \(>k\), after \(\approx 2k\) rows.
• Horizontal tiling because \(k\ll w\) makes \(O(n \lceil k/w\rceil)\) inefficient.

SIMD: Chunking the text

• Split the text in 4 chunks.
• Process 1 chunk per SIMD-lane.
  • The blue blocks are all evaluated at the same time!

Sassy 1: search long text at \({>}1\) Gbp/s

• Params:
  • \(20 \leq |P| \leq 1000\), \(n=10^5\)
  • 1 thread, fwd search only
• Results:
  • k=3: 1.1 Gbp/s
  • k=20: 600 Mbp/s
  • 10\(\times\) Edlib, 128\(\times\) parasail

Applications:

• Crispr off-target searching
  • Many short patterns, fixed reference → Columba (Renders+ 2025)
• Barcode detection for demultiplexing
  • Many short patterns, reads → Batching!

CLI: sassy grep -p <pattern> -k 3 <hg>.fa

Sassy 2: batching \(\gg\) chunking

• Text chunking requires gathering and transposing
• Take pattern suffix of len \(w'\in\{16,32\}\)
  • Random DNA has edit distance \(\approx 45\%\).
  • \(w'>2k+\varepsilon\) ensures few spurious matches.
• One suffix per SIMD lane.
• Post-process if suffix matches with cost \(\leq k\).

Sassy 2: search at 8 Gbp/s

• A: 128 patterns of len 23, 1 thread, \(k\in \{0,3,4\}\). B: \(n=10^5\), \(k=3\).

Sassy 2 Application: Crispr off-target & Barbell

• Search a 23bp guide RNA in a human genome with \(k=3\) in 30ms or amortized >100 Gbp/s on 16 threads.
  • 3.7\(×\) Sassy1, 35\(×\) Edlib
• Scan Nanopore reads for 96 ONT rapid barcodes with \(k=3\) at 1.2 Gbp/s!
  • 4.6\(×\) Sassy1, 45\(×\) Edlib
• Barbell (Beeloo+) is a fast & accurate demultiplexer using Sassy2.
  • Poster tonight!

Outlook: Sassy-GPU

Extras

A*PA: seed heuristic

A*PA: gap-chaining seed heuristic

A*PA: contours

A*PA2: pre-pruning

A*PA2: results (real data)

A*PA2: results (synthetic)