Multiple sequence alignment

Deep dive

How it works

Parts:

• scoring model — score per symbol pair and per gap; gap against gap is 0
• SP score — sum over all row pairs and columns
• distance matrix — e.g. 1 − identity of each pairwise alignment, or k-mer counts
• guide tree — UPGMA (average linkage) or neighbour joining on the distances
• profile — a partial alignment, handled as a sequence of columns
• profile aligner — Needleman–Wunsch where a cell scores one column of A against one column of B

Operations:

1. Compute the N(N−1)/2 pairwise distances.
2. Build the guide tree.
3. From the leaves up, align the two children's profiles. A column pair scores the sum over every row of A × every row of B; a column against a gap column pays the gap cost per residue.
4. Insert the new gap columns into every row of both children; existing gaps never move.
5. The root is the MSA. Optionally refine.

          root = MSA of 4 rows          column score at the root = Σ over 2×2 row pairs
         /          \
   P12 (2 rows)    P34 (2 rows)         P12 = align(s1, s2), P34 = align(s3, s4)
    /    \          /    \
  s1     s2        s3     s4            leaves: single sequences

Variants:

• Iterative refinement (MUSCLE, Edgar 2004): k-mer distances, a progressive pass, then refinement by tree-dependent restricted partitioning (cut a tree edge, re-align the two halves, keep if the score improves).
• Consistency (T-Coffee, Notredame et al. 2000): a library of pairwise alignments, so each merge also weighs how all other sequences align.
• FFT (MAFFT, Katoh et al. 2002): residues become volume and polarity values; an FFT finds homologous regions; FFT-NS-2 is progressive, FFT-NS-i adds refinement.
• HMM–HMM (Clustal Omega, Sievers et al. 2011): mBed guide tree in O(N log N), merges by HHalign on two profile HMMs.

For structure models the MSA comes from a search (profile HMM search): jackhmmer and HHblits (Remmert et al. 2012) collect homologs, and AlphaFold 2 embeds the MSA with pairwise features (Jumper et al. 2021).

Code

Progressive MSA on byte symbols (0 = gap), linear gaps, UPGMA. Tests: letters with indels; word-level log lines whose constant columns give the template.

pub const GAP: u8 = 0; // symbols are non-zero bytes (letters, or ids of words)
pub struct Scoring { pub matched: i32, pub mismatch: i32, pub gap: i32 } // gap < 0, linear

fn sub(a: u8, b: u8, s: &Scoring) -> i32 {
    match (a == GAP, b == GAP) {
        (true, true) => 0,
        (true, false) | (false, true) => s.gap,
        _ => if a == b { s.matched } else { s.mismatch },
    }
}

/// Sum-of-pairs score of one column of A against one column of B (every row against every row).
fn cs(a: &[u8], b: &[u8], s: &Scoring) -> i32 {
    a.iter().map(|&x| b.iter().map(|&y| sub(x, y, s)).sum::<i32>()).sum()
}

fn columns(rows: &[Vec<u8>]) -> Vec<Vec<u8>> {
    (0..rows[0].len()).map(|i| rows.iter().map(|r| r[i]).collect()).collect()
}

/// Needleman-Wunsch on profiles. One row each = plain pairwise alignment. Rows of `a` come first.
pub fn align_profiles(a: &[Vec<u8>], b: &[Vec<u8>], s: &Scoring) -> Vec<Vec<u8>> {
    let (ca, cb) = (columns(a), columns(b));
    let (n, m) = (ca.len(), cb.len());
    let (ga, gb) = (vec![GAP; a.len()], vec![GAP; b.len()]);
    let mut h = vec![vec![0i32; m + 1]; n + 1];
    let mut mv = vec![vec![0u8; m + 1]; n + 1]; // 0 diagonal, 1 gap in B, 2 gap in A
    for i in 0..=n { for j in 0..=m {
        if i + j == 0 { continue; }
        let mut best = (i32::MIN, 0u8);
        if i > 0 && j > 0 { best = (h[i - 1][j - 1] + cs(&ca[i - 1], &cb[j - 1], s), 0); }
        if i > 0 { let v = h[i - 1][j] + cs(&ca[i - 1], &gb, s); if v > best.0 { best = (v, 1); } }
        if j > 0 { let v = h[i][j - 1] + cs(&ga, &cb[j - 1], s); if v > best.0 { best = (v, 2); } }
        (h[i][j], mv[i][j]) = best;
    } }
    let mut out = vec![Vec::new(); a.len() + b.len()];
    let (mut i, mut j) = (n, m);
    while i > 0 || j > 0 {
        let d = mv[i][j];
        let xa = if d != 2 { &ca[i - 1] } else { &ga };
        let xb = if d != 1 { &cb[j - 1] } else { &gb };
        for (r, &c) in xa.iter().chain(xb).enumerate() { out[r].push(c); }
        if d != 2 { i -= 1; }
        if d != 1 { j -= 1; }
    }
    out.iter_mut().for_each(|r| r.reverse());
    out
}

pub fn sp_score(rows: &[Vec<u8>], s: &Scoring) -> i32 {
    let mut t = 0;
    for p in 0..rows.len() { for q in p + 1..rows.len() {
        t += rows[p].iter().zip(&rows[q]).map(|(&x, &y)| sub(x, y, s)).sum::<i32>();
    } }
    t
}

/// Progressive MSA: pairwise distances, UPGMA merge order, profile merges up the tree.
pub fn progressive(seqs: &[Vec<u8>], s: &Scoring) -> Vec<Vec<u8>> {
    let k = seqs.len();
    let mut d = vec![vec![0.0f64; k]; k];
    for p in 0..k { for q in p + 1..k {
        let al = align_profiles(&seqs[p..=p], &seqs[q..=q], s);
        let same = al[0].iter().zip(&al[1]).filter(|(x, y)| x == y).count();
        d[p][q] = 1.0 - same as f64 / al[0].len() as f64; // 1 - identity
        d[q][p] = d[p][q];
    } }
    let mut cl: Vec<(Vec<usize>, Vec<Vec<u8>>)> = (0..k).map(|i| (vec![i], vec![seqs[i].clone()])).collect();
    while cl.len() > 1 {
        let (mut best, mut bp, mut bq) = (f64::MAX, 0, 1);
        for p in 0..cl.len() { for q in p + 1..cl.len() {
            let sum: f64 = cl[p].0.iter().map(|&x| cl[q].0.iter().map(|&y| d[x][y]).sum::<f64>()).sum();
            let avg = sum / (cl[p].0.len() * cl[q].0.len()) as f64; // UPGMA: average linkage
            if avg < best { (best, bp, bq) = (avg, p, q); }
        } }
        let (mb, rb) = cl.remove(bq);
        let (ma, ra) = cl.remove(bp);
        cl.push(([ma, mb].concat(), align_profiles(&ra, &rb, s)));
    }
    let (members, rows) = cl.pop().unwrap();
    let mut out = vec![Vec::new(); k];
    for (r, i) in rows.into_iter().zip(members) { out[i] = r; }
    out
}

#[cfg(test)]
mod tests {
    use super::*;
    const S: Scoring = Scoring { matched: 2, mismatch: -1, gap: -2 };

    #[test]
    fn letters_with_indels() {
        let seqs: Vec<Vec<u8>> = ["GATTACA", "GATACA", "GATTTACA", "GCTTACA"].iter().map(|x| x.bytes().collect()).collect();
        let al = progressive(&seqs, &S);
        assert!(al.iter().all(|r| r.len() == 8)); // 8 = longest input: all indels share one column
        for (r, q) in al.iter().zip(&seqs) { assert_eq!(&r.iter().copied().filter(|&c| c != GAP).collect::<Vec<_>>(), q); }
        assert_eq!(sp_score(&al, &S), 57); // measured
    }

    #[test]
    fn log_line_template() {
        let lines = ["user alice logged in from host a", "user bob logged in from host b",
                     "user carol logged in via ssh from host c", "user dave logged in from host d"];
        let mut vocab: Vec<&str> = vec![];
        let toks: Vec<Vec<u8>> = lines.iter().map(|l| l.split(' ').map(|w| {
            if !vocab.contains(&w) { vocab.push(w); }
            vocab.iter().position(|&v| v == w).unwrap() as u8 + 1
        }).collect()).collect();
        let al = progressive(&toks, &S);
        let t: Vec<&str> = (0..al[0].len())
            .map(|c| if al.iter().all(|r| r[c] == al[0][c]) { vocab[al[0][c] as usize - 1] } else { "*" })
            .collect();
        assert_eq!(t.join(" "), "user * logged in * * from host *");
    }
}

cargo test --release passes both; the letters case gives 8 columns and SP 57 (padded without alignment: 15).

Complexity and constants

• Distances: N(N−1)/2 alignments of O(L²).
• Profile merge: O(L²·|A|·|B|). Over the tree Σ|A||B| = N(N−1)/2, so merging is also O(N²L²).
• This UPGMA is O(N³); mBed builds trees in O(N log N) (Sievers et al. 2011).
• Exact SP: O(L^N) cells; Carrillo and Lipman (1988) prune the lattice using pairwise alignments.
• Memory: one (L+1)² i32 table per merge; L = 10⁴ is order of 400 MB.

Efficiency

• Cost: the code above, one thread, on copies of a random DNA root with 5 % deletions, 5 % insertions and 5 % substitutions per position, measured on this Mac (Apple M4 Pro): N = 10, L = 200 takes 13 ms; N = 20, 52 ms; N = 40, 205 ms; N = 20, L = 400 takes 205 ms. Doubling N or L costs about 4×, as O(N²L²) says. Bound: DP cells; memory for long L.
• Levers:
  • k-mer counts instead of DP for distances — MUSCLE aligns 5000 sequences of length 350 in 7 min on a 2004 desktop (Edgar 2004) — large N
  • mBed guide tree — O(N log N) (Sievers et al. 2011) — N in the thousands
  • SIMD bit-parallel Myers distances, bisecting k-means trees — Kalign 3 (Lassmann 2020) — very large sets
  • FFT to find homologous stretches — FFT-NS-i over 100 times faster than T-Coffee above 60 sequences (Katoh et al. 2002) — long, similar sequences
• Naive → standard → state of the art:
  • exact N-dimensional DP · O(L^N)
  • pairwise distances, UPGMA, profile merges · O(N²L²)
  • k-mer or embedding trees, SIMD, HMM profiles · near-linear guide tree, merges dominate

Implementation notes

• Large alphabets: keep per-column residue counts and a precomputed matrix product; a cell costs O(σ), not |A||B|.
• Real data wants affine, position-dependent gap penalties; the code uses linear gaps.
• Weight sequences so near-duplicates do not dominate SP (Clustal W, Thompson et al. 1994).
• Symbol 0 is the gap; map words to ids from 1; above 255 symbols use u16.
• Alignments are not unique: test properties (rows ungap to inputs, equal length, SP), not exact strings.

Pitfalls

• Greedy merges freeze mistakes: a misplaced gap in an early pair never moves.
• Unrelated or fragmentary inputs give meaningless columns.
• On real logs, normalise tokens first or numbers and ids all become slots.

Variants and alternatives

In production

• AlphaFold 2's pipeline: jackhmmer and HHblits MSAs as network input (Jumper et al. 2021).
• MAFFT, MUSCLE, Clustal Omega, Kalign for gene and protein families; CollateX for manuscript collation.

Implementations

Sources

• Feng, Doolittle, Progressive sequence alignment as a prerequisite to correct phylogenetic trees, J Mol Evol 25, 1987. paper
• Carrillo, Lipman, The multiple sequence alignment problem in biology, SIAM J Appl Math 48, 1988. paper
• Wang, Jiang, On the complexity of multiple sequence alignment, J Comput Biol 1, 1994. paper
• Thompson, Higgins, Gibson, CLUSTAL W, Nucleic Acids Res 22, 1994. paper
• Notredame, Higgins, Heringa, T-Coffee, J Mol Biol 302, 2000. paper
• Katoh et al., MAFFT, Nucleic Acids Res 30, 2002. paper
• Edgar, MUSCLE, Nucleic Acids Res 32, 2004. paper
• Sievers et al., Clustal Omega, Mol Syst Biol 7, 2011. paper
• Remmert et al., HHblits, Nat Methods 9, 2012. paper
• Lassmann, Kalign 3, Bioinformatics 36, 2020. paper · code
• Jumper et al., Highly accurate protein structure prediction with AlphaFold, Nature 596, 2021. paper
• Edgar, Muscle5, Nat Commun 13, 2022. paper
• Haentjens Dekker et al., Computer-supported collation of modern manuscripts, Digital Scholarship in the Humanities 30, 2015. paper · code
• Fiscus, A post-processing system to yield reduced word error rates (ROVER), IEEE ASRU Workshop 1997. paper