A history of pairwise alignment

My blog: curiouscoding.nl Link to heading

Pairwise alignment:
finding differences between strings Link to heading

Covid – \(\alpha\), December 2020 Link to heading

Covid – \(\omicron\), December 2021 Link to heading

Pairwise alignment Link to heading

• Find the mutations between two sequences

Dynamic programming Link to heading

Dynamic programming Link to heading

Dynamic programming Link to heading

Dynamic programming Link to heading

Dynamic programming Link to heading

50 000 TB of sequenced DNA → it must be fast Link to heading

Goal:
Fast code Link to heading

Goal:
High throughput code Link to heading

Goal:
Optimal throughput code Link to heading

What is high troughput code? Link to heading

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! Link to heading

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

Can we do better than quadratic? Link to heading

No! (not much anyway) Link to heading

• 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 Link to heading

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

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

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

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

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

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

A*PA heuristics Link to heading

A*PA – not quite linear, and terrible efficiency Link to heading

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

Band Doubling – \(O(ns)\) Link to heading

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

Optimizations Link to heading

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 Link to heading

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

A*PA: comparison Link to heading

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

Alignment modes Link to heading

Semi-global variants Link to heading

… Sassy: 100bp semi-global alignment Link to heading

Approximate String Matching Link to heading

• 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)\) Link to heading

• 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 Link to heading

• 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 Link to heading

• 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:

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

Sassy 2: batching \(\gg\) chunking Link to heading

• 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 Link to heading

• 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 Link to heading

• 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 Link to heading

Extras Link to heading

A*PA: seed heuristic Link to heading

A*PA: gap-chaining seed heuristic Link to heading

A*PA: contours Link to heading

A*PA2: pre-pruning Link to heading

A*PA2: results (real data) Link to heading

A*PA2: results (synthetic) Link to heading