Smith-Waterman local alignment

Deep dive

How it works

Parts:

• scoring model — s(x,y) match/mismatch or a substitution table (BLOSUM62 for proteins; +match/-mismatch for DNA), gap open ≥ extend > 0
• H[i][j] — best score of an alignment ending at a[i-1], b[j-1], floored at 0
• E[i][j] — best score ending in a gap that consumes b[j-1] (horizontal run)
• F[i][j] — best score ending in a gap that consumes a[i-1] (vertical run)
• best cell — running maximum of H, start of the traceback

Operations:

1. Init H = 0 on row 0 and column 0; E = F = -inf.
2. For each cell: E = max(H[i][j-1] - open, E[i][j-1] - extend); F = max(H[i-1][j] - open, F[i-1][j] - extend).
3. H = max(0, H[i-1][j-1] + s(a[i-1], b[j-1]), E, F). The 0 is the whole difference from global alignment: a bad prefix is dropped.
4. Track the maximum H (and its cell) over the table.
5. Traceback from that cell in state H: stop when H == 0; if H came from the diagonal emit a pair; else switch to E or F and emit gaps while the cell was reached by extending (E == E_prev - extend); leave the state where the gap was opened from H.
6. Reverse the emitted pairs. Start = the cell where the walk stopped.

 H[i-1][j-1] ──s(a_i,b_j)──┐
 H[i][j-1] ─open─► E ──────┼─► H[i][j] = max(0, diag, E, F)
 H[i-1][j] ─open─► F ──────┘   (E, F also extend themselves)
 best = max H anywhere ─► walk back through H/E/F until H = 0

Statistics (from memory, unverified): for ungapped local scores with negative expected pair score, chance hits follow E = K·m·n·e^(−λS) (Karlin and Altschul 1990); λ, K depend on the scoring system and letter frequencies and are fitted for gapped scoring.

Code

Generic over the symbol type, so bytes and tokens use the same function. Gap of length k costs open + (k-1)·extend (EMBOSS and WU-BLAST convention).

pub struct Scoring { pub matc: i32, pub mismatch: i32, pub open: i32, pub extend: i32 }
pub struct Alignment {
    pub score: i32,
    pub start: (usize, usize), // first aligned cell (0-based) in a, b
    pub end: (usize, usize),   // one past the last aligned cell
    pub pairs: Vec<(Option<usize>, Option<usize>)>, // None = gap
}
const NEG: i32 = i32::MIN / 2; // "minus infinity" that survives a subtraction

pub fn smith_waterman<T: PartialEq>(a: &[T], b: &[T], sc: &Scoring) -> Alignment {
    let (n, m) = (a.len(), b.len());
    let w = m + 1;
    let (mut h, mut e, mut f) = (vec![0i32; (n + 1) * w], vec![NEG; (n + 1) * w], vec![NEG; (n + 1) * w]);
    let sub = |x: &T, y: &T| if x == y { sc.matc } else { sc.mismatch };
    let (mut best, mut bi, mut bj) = (0, 0, 0);
    for i in 1..=n {
        for j in 1..=m {
            let k = i * w + j;
            e[k] = (h[k - 1] - sc.open).max(e[k - 1] - sc.extend); // consumes b[j-1]
            f[k] = (h[k - w] - sc.open).max(f[k - w] - sc.extend); // consumes a[i-1]
            h[k] = 0.max(h[k - w - 1] + sub(&a[i - 1], &b[j - 1])).max(e[k]).max(f[k]);
            if h[k] > best { (best, bi, bj) = (h[k], i, j); }
        }
    }
    // Traceback from the maximum cell until H hits 0. state 0=H, 1=E, 2=F.
    let (mut i, mut j, mut st) = (bi, bj, 0);
    let mut pairs = Vec::new();
    loop {
        let k = i * w + j;
        match st {
            0 if h[k] == 0 => break,
            0 if h[k] == h[k - w - 1] + sub(&a[i - 1], &b[j - 1]) => {
                pairs.push((Some(i - 1), Some(j - 1)));
                (i, j) = (i - 1, j - 1);
            }
            0 => st = if h[k] == e[k] { 1 } else { 2 },
            1 => { // leave E when the gap was opened here (came from H)
                pairs.push((None, Some(j - 1)));
                if e[k] != e[k - 1] - sc.extend { st = 0; }
                j -= 1;
            }
            _ => {
                pairs.push((Some(i - 1), None));
                if f[k] != f[k - w] - sc.extend { st = 0; }
                i -= 1;
            }
        }
    }
    pairs.reverse();
    Alignment { score: best, start: (i, j), end: (bi, bj), pairs }
}

/// Score only, O(len(b)) memory: keep one row of H and F, E is a running scalar.
pub fn score_only<T: PartialEq>(a: &[T], b: &[T], sc: &Scoring) -> i32 {
    let m = b.len();
    let (mut h, mut f, mut best) = (vec![0i32; m + 1], vec![NEG; m + 1], 0);
    for x in a {
        let (mut diag, mut left, mut e) = (0, 0, NEG);
        for j in 1..=m {
            e = (left - sc.open).max(e - sc.extend);
            f[j] = (h[j] - sc.open).max(f[j] - sc.extend);
            let v = 0.max(diag + if *x == b[j - 1] { sc.matc } else { sc.mismatch }).max(e).max(f[j]);
            (diag, h[j], left) = (h[j], v, v);
            best = best.max(v);
        }
    }
    best
}

Compiled and run on this Mac with these asserts passing: hand-checked TGTTACGG vs GGTTGACTA (match 3, mismatch -3, gap 2: score 13, GTT-AC over GTTGAC); one 5-long gap scoring 32 - (4+4) = 24; identical strings; AAAA vs CCCC (score 0, no pairs); a 12-base core inside unrelated flanks (core only); score_only equal to the full score; and two code-token streams where renamed variables become mismatches (score 25, 17 pairs).

Complexity and constants

• Time O(nm); space 12·nm bytes for the 3 tables with traceback, 8·min(n,m) bytes (+ a few scalars) for score only.
• Measured on this Mac (Apple M4 Pro), --release, one thread, 1,000 × 10,000 random DNA, match 2, mismatch -3, open 5, extend 2, 3 runs each: score_only 0.39-0.55 GCUPS (18-25 ms); full tables + traceback 0.31-0.61 GCUPS (16-33 ms). Noisy at this size: order of 0.5 GCUPS.
• Full tables at 1,000 × 10,000 hold 3 × 10 M × 4 B = 120 MB; at 100 k × 100 k, 120 GB (arithmetic), hence score-only and linear-space variants.
• Each cell needs left, up and diagonal, so the row loop is a serial chain through E and H; SIMD layouts exist to break it.

Efficiency

• Cost: one cell update, measured on this Mac (Apple M4 Pro), scalar: about 2 ns (0.5 GCUPS, run above). Bound: compute and the serial E chain; memory only for full tables.
• Levers:
  • striped SIMD — Farrar 2007 abstract: 2-8x over other SIMD implementations, >3 GCUPS on a 2.0 GHz Xeon — score-only passes over many sequences.
  • banding / X-drop — cells O(nm) to O(n·band) (algorithmic) — endpoints known near a diagonal, as after a seed hit.
  • seed filter first — DP only on candidate windows — long references; see Minimizer k-mer index.
  • score only first, traceback for survivors — drops the 12 B/cell tables — most candidates rejected.
  • 8/16-bit saturating lanes with overflow fallback — more lanes per register — short scores in SIMD code.
• Naive → standard → state of the art:
  • full tables, traceback for every pair · order of 0.5 GCUPS here, 12 B/cell
  • rolling rows, traceback on winners only · same speed, O(min) memory
  • striped or inter-sequence SIMD, banded/adaptive for long reads, GPU batches · figures from the cited papers, not measured here

Implementation notes

• Global (Needleman-Wunsch): no floor at 0, start at the corner. Semi-global/overlap: free end gaps, no floor. Local: floor at 0. Same loop; init, floor and where the maximum is taken differ.
• Gap conventions differ: first gap cell costs open (EMBOSS, WU-BLAST) or open + extend (NCBI-BLAST), per the EMBOSS water manual. Fix one in your test fixtures.
• Linear-space traceback: Hirschberg 1975 splits at the middle row; Myers and Miller 1988 give the affine (Gotoh) form.
• Striped SIMD (Farrar): query column split into segLen vectors so the vertical dependency is a vector op, plus a short corrective loop. Inter-sequence SIMD puts one database sequence per lane.
• Tie-breaking decides among equal-score alignments; keep it deterministic.

Pitfalls

• Non-negative expected pair score: local alignments degenerate to near-global.
• E-values with λ, K not fitted for the matrix and gap costs in use.
• NEG = i32::MIN then subtracting overflows; use MIN/2 or saturating math.
• Comparing raw scores across different query lengths; normalise or use E-values.
• Low-complexity runs (poly-A, repeated tokens) score high; mask them first.

Variants and alternatives

In production

• BLAST (Altschul et al. 1990) searches by seed-then-extend rather than full DP against every database sequence.
• Read mapping: minimap2 includes the KSW2 aligner for base-level alignment (ksw2 README).
• EMBOSS water; parasail (diagonal, blocked, striped and prefix-scan SIMD variants).

Implementations

Sources

• Temple F. Smith, Michael S. Waterman, Identification of common molecular subsequences, J Mol Biol 147(1):195-197, 1981. paper
• Osamu Gotoh, An improved algorithm for matching biological sequences, J Mol Biol 162(3):705-708, 1982. paper
• Saul B. Needleman, Christian D. Wunsch, A general method applicable to the search for similarities in the amino acid sequence of two proteins, J Mol Biol 48(3):443-453, 1970. paper
• Michael Farrar, Striped Smith-Waterman speeds database searches six times over other SIMD implementations, Bioinformatics 23(2):156-161, 2007 (doi 10.1093/bioinformatics/btl582; publisher blocks scripted access, record confirmed via Crossref).
• Stephen F. Altschul, Warren Gish, Webb Miller, Eugene W. Myers, David J. Lipman, Basic local alignment search tool, J Mol Biol 215(3):403-410, 1990. paper
• Heng Li, Minimap2: pairwise alignment for nucleotide sequences, Bioinformatics 34(18):3094-3100, 2018 (doi 10.1093/bioinformatics/bty191). code · ksw2
• Samuel Karlin, Stephen F. Altschul, Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes, PNAS 87(6):2264-2268, 1990 (doi 10.1073/pnas.87.6.2264).
• Steven Henikoff, Jorja G. Henikoff, Amino acid substitution matrices from protein blocks, PNAS 89(22):10915-10919, 1992 (doi 10.1073/pnas.89.22.10915; BLOSUM).
• Daniel S. Hirschberg, A linear space algorithm for computing maximal common subsequences, CACM 18(6):341-343, 1975 (doi 10.1145/360825.360861). Eugene W. Myers, Webb Miller, Optimal alignments in linear space, CABIOS 4(1):11-17, 1988 (doi 10.1093/bioinformatics/4.1.11).
• Jeff Daily, Parasail: SIMD C library for global, semi-global, and local pairwise sequence alignments, BMC Bioinformatics 17:81, 2016 (doi 10.1186/s12859-016-0930-z). code
• EMBOSS water manual (gap conventions). docs