Computational Biology Algorithms, in depth: local alignment, minimizers, seed-chain-extend mapping, HMMs and neighbor joining

The map: from reads to answers

The figure follows a resequencing analysis: reads from a sample are mapped to a reference genome, and differences become variant calls. Mapping is seed, chain, extend: find short exact matches using an index, keep groups of matches that lie along a consistent diagonal, then run a full alignment only in the small region those matches point to. When there is no reference, reads are assembled instead, usually through a de Bruijn or overlap graph. Downstream, probabilistic models annotate the sequence (genes, CpG islands, protein families) and phylogenetic methods relate sequences across species.

The recurring idea is that the exact algorithm is quadratic, and the data makes quadratic impossible: aligning a 150-base read against 3.1 billion positions by full dynamic programming is around 4.7 × 10¹¹ cell updates, for one read. Every practical method spends most of its effort deciding where not to look.

Local alignment: Smith-Waterman with affine gaps

Global alignment (Needleman-Wunsch) aligns two sequences end to end. Biology more often asks a local question: which part of this read, or this protein, matches something in that one? Smith-Waterman answers it with one change to the recurrence: every cell may also take the value 0, meaning "start a new alignment here", and the answer is the maximum cell anywhere in the table rather than the bottom-right corner.

One five-base insertion is far likelier than five separate ones, so affine scoring charges an opening penalty plus a smaller extension penalty. Gotoh's formulation keeps it O(n·m) with three matrices: ending in a match or mismatch, or inside a gap in either sequence.

NEG = float("-inf")

def smith_waterman_affine(a, b, match=2, mismatch=-3, gap_open=-5, gap_extend=-2):
    """Gotoh local alignment. A gap of length L costs gap_open + (L-1)*gap_extend."""
    n, m = len(a), len(b)
    H = [[0] * (m + 1) for _ in range(n + 1)]       # best local score ending at (i, j)
    E = [[NEG] * (m + 1) for _ in range(n + 1)]     # ... ending in a gap in a
    F = [[NEG] * (m + 1) for _ in range(n + 1)]     # ... ending in a gap in b
    best, where = 0, (0, 0)
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            E[i][j] = max(E[i][j - 1] + gap_extend, H[i][j - 1] + gap_open)
            F[i][j] = max(F[i - 1][j] + gap_extend, H[i - 1][j] + gap_open)
            s = match if a[i - 1] == b[j - 1] else mismatch
            H[i][j] = max(0, H[i - 1][j - 1] + s, E[i][j], F[i][j])
            if H[i][j] > best:
                best, where = H[i][j], (i, j)
    return best, where

print(smith_waterman_affine("TTGACACCCTCCCAATT", "ACCTCCTAAT"))   # (13, (16, 10))

The best local alignment scores 13 and ends at position 16 of the first sequence and 10 of the second. It is CCTCCCAAT against CCTCCTAAT: eight matches and one mismatch, 16 − 3 = 13. A gapped alternative that starts earlier, ACCCTCCCAAT against ACC-TCCTAAT, scores only 10: opening the gap costs 5 and buys just one extra match, worth 2. Change the penalties and the winner changes, which is why scoring parameters, usually derived from substitution matrices such as BLOSUM for proteins, matter as much as the algorithm. For the global version, traceback and banding in depth, see Needleman-Wunsch and edit distance.

Minimizers: sampling k-mers without losing matches

Indexing every k-mer of a genome is large, and most of them are redundant because neighbouring k-mers overlap. A minimizer scheme slides a window of w consecutive k-mers along the sequence and keeps only the smallest one in each window, under some fixed order. Two properties make this useful. Any two sequences that share a stretch of at least w + k − 1 bases are guaranteed to share a minimizer, so sampling never loses a long enough match. And the choice depends only on local content, so the same substring in the read and the reference picks the same k-mers.

def minimizers(seq, k, w):
    """(position, k-mer) of the smallest k-mer in each window of w consecutive k-mers.
    Real tools order k-mers by a hash and use the canonical strand; plain string
    order keeps this example checkable by hand."""
    kmers = [seq[i:i + k] for i in range(len(seq) - k + 1)]
    out = []
    for s in range(len(kmers) - w + 1):
        pos = min(range(s, s + w), key=lambda i: (kmers[i], i))
        if not out or out[-1][0] != pos:
            out.append((pos, kmers[pos]))
    return out

print(minimizers("GATTACAGATTACA", 3, 4))
# [(1, 'ATT'), (4, 'ACA'), (6, 'AGA'), (8, 'ATT'), (11, 'ACA')]

The sequence has 12 three-mers, but only 5 are kept, and the repeated GATTACA selects ATT and ACA both times. With random-looking input and a hash order, the expected density is about 2/(w+1), so large windows give large savings. Lexicographic order, as used here for readability, is a poor choice in practice: it over-samples low-complexity k-mers such as runs of A. The same sampling idea, with a different selection rule, powers MinHash sketches that estimate the similarity of whole genomes from a few thousand hashes.

Read mapping: seed, chain, extend

Seeds are short exact matches between a read and the reference, found either by looking up the read's minimizers in a hash table of the reference's minimizers, or by backward search in an FM-index built on the Burrows-Wheeler transform, which finds the range (and count) of all occurrences in time proportional to the pattern length; reporting each position costs extra via a sampled suffix array. The Burrows-Wheeler transform and suffix array construction pages cover those indexes.

A read hits many places by chance, especially in repeats. Chaining filters them: a true alignment produces anchors that increase together in both the reference and the read, along roughly the same diagonal. Dynamic programming over the anchors finds the best such chain.

def chain(anchors, max_gap=50):
    """anchors: (ref_pos, read_pos, length), sorted by ref_pos.
    Highest-scoring co-linear chain by O(n^2) dynamic programming."""
    n = len(anchors)
    score = [a[2] for a in anchors]
    prev = [-1] * n
    for i in range(n):
        ri, qi, li = anchors[i]
        for j in range(i):
            rj, qj, lj = anchors[j]
            dr, dq = ri - rj, qi - qj
            if dr <= 0 or dq <= 0 or max(dr, dq) > max_gap:
                continue
            gain = min(li, dr, dq) - abs(dr - dq)    # new bases minus diagonal drift
            if score[j] + gain > score[i]:
                score[i], prev[i] = score[j] + gain, j
    end = max(range(n), key=score.__getitem__)
    path = []
    while end != -1:
        path.append(anchors[end]); end = prev[end]
    return path[::-1], max(score)

Given anchors at reference/read positions (100, 0), (118, 17), (140, 39), (160, 60) and a stray hit at (300, 20), all of length 15, the chain keeps the first four with a score of 58 and drops the stray one, which is off the diagonal. Production mappers bound how many predecessors each anchor examines, making this near-linear. Finally, a banded Smith-Waterman extends around the chain, filling only cells near its diagonal.

Hidden Markov models and Viterbi

Many biological questions are segmentation questions: which parts of this sequence are genes, which are CpG islands (regions rich in C and G that often mark gene promoters), which residues of a protein belong to a known family? A hidden Markov model describes a sequence as generated by a walk through hidden states, each with its own letter frequencies, and the Viterbi algorithm finds the most likely walk for an observed sequence by dynamic programming over positions and states, in O(length × states²) time.

import math

def viterbi(obs, states, start, trans, emit):
    """Most likely hidden-state path, in log space to avoid underflow."""
    V = [{s: math.log(start[s]) + math.log(emit[s][obs[0]]) for s in states}]
    back = [{}]
    for t in range(1, len(obs)):
        V.append({}); back.append({})
        for s in states:
            prev = max(states, key=lambda r: V[t - 1][r] + math.log(trans[r][s]))
            V[t][s] = V[t - 1][prev] + math.log(trans[prev][s]) + math.log(emit[s][obs[t]])
            back[t][s] = prev
    last = max(states, key=lambda s: V[-1][s])
    path = [last]
    for t in range(len(obs) - 1, 0, -1):
        path.append(back[t][path[-1]])
    return path[::-1], V[-1][last]

states = ("CpG", "bg")
start = {"CpG": 0.5, "bg": 0.5}
trans = {"CpG": {"CpG": 0.9, "bg": 0.1}, "bg": {"CpG": 0.1, "bg": 0.9}}
emit = {"CpG": {"A": .15, "C": .35, "G": .35, "T": .15},
        "bg":  {"A": .30, "C": .20, "G": .20, "T": .30}}
# ATATTACGCGGCGCCGATAATTA  ->  ......CCCCCCCCCC.......

This toy model has two states: a C/G-rich island state and a background state, each sticky with probability 0.9. On ATATTACGCGGCGCCGATAATTA Viterbi labels positions 6 to 15, exactly the CGCGGCGCCG stretch, as island. Work in log space; multiplying hundreds of probabilities below 1 underflows to zero in floating point. Real CpG models track dinucleotides, since the signal is the frequency of the pair CG, not of C and G separately, and gene finders and profile HMMs use many more states, but the decoding algorithm is exactly this one.

Phylogenetic trees: neighbor joining

Once sequences are aligned, their pairwise distances suggest how they are related. Neighbor joining builds a tree from a distance matrix without assuming every lineage evolves at the same rate. At each step it joins the pair that minimises a corrected score Q(i, j) = (n − 2)·d(i, j) − r(i) − r(j), where r is a row sum; the correction stops it from merely pairing the two closest taxa, which goes wrong when one lineage evolves fast. It computes branch lengths for the joined pair, replaces them with a new node, and repeats. Each step is O(n²), so the whole tree is O(n³).

def neighbor_joining(names, D):
    """Saitou-Nei neighbor joining on a symmetric distance matrix. Returns Newick."""
    nodes, D = list(names), [row[:] for row in D]
    while len(nodes) > 2:
        n = len(nodes)
        r = [sum(row) for row in D]
        q, i, j = min(((n - 2) * D[i][j] - r[i] - r[j], i, j)
                      for i in range(n) for j in range(i + 1, n))
        li = 0.5 * D[i][j] + (r[i] - r[j]) / (2 * (n - 2))
        lj = D[i][j] - li
        joined = f"({nodes[i]}:{li:g},{nodes[j]}:{lj:g})"
        du = [0.5 * (D[i][k] + D[j][k] - D[i][j]) for k in range(n)]
        keep = [k for k in range(n) if k not in (i, j)]
        D = [[D[a][b] for b in keep] + [du[a]] for a in keep] + [[du[a] for a in keep] + [0]]
        nodes = [nodes[k] for k in keep] + [joined]
    return f"({nodes[0]},{nodes[1]}:{D[0][1]:g});"

On the standard five-taxon example (distances a–b 5, a–c 9, a–d 9, a–e 8, b–c 10, b–d 10, b–e 9, c–d 8, c–e 7, d–e 3) the run returns ((c:4,(a:2,b:3):3),(d:2,e:1):2);: a and b join first with branch lengths 2 and 3, d and e pair with lengths 2 and 1, and c hangs off the a–b side. Because the matrix is exactly additive, the tree reproduces every input distance. Real distances are noisy, so check trees with bootstrap support.

Operational guidance

• Know your data's error model. Short reads mostly carry substitutions; long reads carry far more indels, which changes seed length, gap penalties and chaining tolerance.
• Mask or down-weight repeats. A seed in a repeat hits thousands of places; mappers cap or skip high-frequency seeds and report a mapping quality that tells you how unique the placement was. Filter on it.
• Remember both strands. DNA is double-stranded, so a read may come from the reverse complement. Index canonical k-mers (the smaller of a k-mer and its reverse complement) or search both orientations.
• Keep exact methods as oracles. Run full Smith-Waterman on a sample of reads to measure what your heuristic misses, and record the reference build, index parameters and scoring scheme with every result.

Failure modes

Trade-offs

Exact dynamic programming is optimal for its scoring model and quadratic in cost; seed-chain-extend is fast but can miss divergent matches whose seeds never hit. Smaller k and denser sampling raise sensitivity and cost. FM-indexes are compact but awkward to update; minimizer hash tables are simple to build but larger. HMMs are only as good as their parameters, and neighbor joining trades accuracy for speed against likelihood methods. For the no-reference branch of the map, read de Bruijn graph assembly.

What to do next

1. Implement Smith-Waterman with affine gaps, add a traceback, and confirm the 13-point alignment in the example.
2. Add a band of width b around a given diagonal and measure how many cells it saves on a 1,000-base pair.
3. Write the minimizer function with a hash order and canonical k-mers, and measure its density against 2/(w+1).
4. Build a toy mapper: index a 100 kb reference by minimizers, chain read hits, extend with your banded aligner, and check placements on simulated reads.
5. Turn the CpG HMM into a dinucleotide model and decode a real promoter sequence.
6. Rebuild the five-taxon neighbor-joining tree by hand, then on your own distance matrix.