smith-waterman-affine
The implementation is in ~/smith_waterman_affine.py. It implements the full Smith-Waterman local alignment algorithm with an affine gap penalty model using three DP matrices:
Three matrices:
- M[i][j] — score ending with residue seq1[i-1] aligned to seq2[j-1]
- Ix[i][j] — score ending with a gap in seq1 (i.e., seq2[j-1] aligned to a gap)
- Iy[i][j] — score ending with a gap in seq2 (i.e., seq1[i-1] aligned to a gap)
Recurrence relations (all take max(0, ...) for local alignment):
M[i][j] = max(0, M[i-1][j-1] + S(a_i,b_j), Ix[i-1][j-1] + S(a_i,b_j), Iy[i-1][j-1] + S(a_i,b_j))
Ix[i][j] = max(0, M[i][j-1] + gap_open, Ix[i][j-1] + gap_extend)
Iy[i][j] = max(0, M[i-1][j] + gap_open, Iy[i-1][j] + gap_extend)
Traceback: follows the maximum-score path back from best_i, best_j by checking which predecessor cell (with the same matrix type) yielded the current value. A safety bound of n + m steps prevents infinite loops.
BLOSUM62: a full 24×24 matrix (20 standard amino acids + B, Z, X, *) is included as a lookup table.
Key code excerpt (core DP loop):
for i in range(1, n + 1):
for j in range(1, m + 1):
s = subst_fn(seq1[i - 1], seq2[j - 1])
M[i, j] = max(0.0,
M[i - 1, j - 1] + s,
Ix[i - 1, j - 1] + s,
Iy[i - 1, j - 1] + s)
Ix[i, j] = max(0.0,
M[i, j - 1] + gap_open,
Ix[i, j - 1] + gap_extend)
Iy[i, j] = max(0.0,
M[i - 1, j] + gap_open,
Iy[i - 1, j] + gap_extend)
val = max(M[i, j], Ix[i, j], Iy[i, j])
if val > best_score:
best_score = val
best_i, best_j = i, j
All 24 test assertions pass across 12 test groups: | Test Group | Assertions | Status | |---|---|---| | Identical sequences | 1 | ✓ | | Empty sequences | 3 | ✓ | | No match possible | 2 | ✓ | | Gap penalty behavior | 2 | ✓ | | BLOSUM62 scores | 3 | ✓ | | Traceback correctness | 3 | ✓ | | Local alignment subset | 1 | ✓ | | Multiple alignments | 2 | ✓ | | Single character sequences | 2 | ✓ | | Reversibility (score symmetric) | 1 | ✓ | | Affine vs linear gap equivalence | 1 | ✓ | | Real protein example | 3 | ✓ | **Edge cases tested:** - **Empty sequences** (seq1 empty, seq2 empty, both empty) → score 0 - **No match possible** (AAAA vs CCCC with very high gap penalties) → score 0 - **Single characters** (A vs A = match, A vs B = mismatch) - **Affine gap penalties** (opening vs extension costs behave correctly) - **Reversibility** (swapping sequences yields identical score) - **Non-standard AA** (`'O'` raises `KeyError` as expected — not in BLOSUM62) - **`X` (any amino acid)** works correctly (X vs X = −1 in BLOSUM62) **Example output** from a canonical globin alignment: ``` Score: 90.0 Aln1: LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF Aln2: LTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRFFESF ``` ---
{"model": "pi", "problem_class": "smith-waterman-affine", "result": "passed", "tests": 24}The implementation is in ~/smith_waterman_affine.py. It implements the full Smith-Waterman local alignment algorithm with an affine gap penalty model using three DP matrices:
Three matrices:
- M[i][j] — score ending with residue seq1[i-1] aligned to seq2[j-1]
- Ix[i][j] — score ending with a gap in seq1 (i.e., seq2[j-1] aligned to a gap)
- Iy[i][j] — score ending with a gap in seq2 (i.e., seq1[i-1] aligned to a gap)
Recurrence relations (all take max(0, ...) for local alignment):
M[i][j] = max(0, M[i-1][j-1] + S(a_i,b_j), Ix[i-1][j-1] + S(a_i,b_j), Iy[i-1][j-1] + S(a_i,b_j))
Ix[i][j] = max(0, M[i][j-1] + gap_open, Ix[i][j-1] + gap_extend)
Iy[i][j] = max(0, M[i-1][j] + gap_open, Iy[i-1][j] + gap_extend)
Traceback: follows the maximum-score path back from best_i, best_j by checking which predecessor cell (with the same matrix type) yielded the current value. A safety bound of n + m steps prevents infinite loops.
BLOSUM62: a full 24×24 matrix (20 standard amino acids + B, Z, X, *) is included as a lookup table.
Key code excerpt (core DP loop):
for i in range(1, n + 1):
for j in range(1, m + 1):
s = subst_fn(seq1[i - 1], seq2[j - 1])
M[i, j] = max(0.0,
M[i - 1, j - 1] + s,
Ix[i - 1, j - 1] + s,
Iy[i - 1, j - 1] + s)
Ix[i, j] = max(0.0,
M[i, j - 1] + gap_open,
Ix[i, j - 1] + gap_extend)
Iy[i, j] = max(0.0,
M[i - 1, j] + gap_open,
Iy[i - 1, j] + gap_extend)
val = max(M[i, j], Ix[i, j], Iy[i, j])
if val > best_score:
best_score = val
best_i, best_j = i, j
All 24 test assertions pass across 12 test groups: | Test Group | Assertions | Status | |---|---|---| | Identical sequences | 1 | ✓ | | Empty sequences | 3 | ✓ | | No match possible | 2 | ✓ | | Gap penalty behavior | 2 | ✓ | | BLOSUM62 scores | 3 | ✓ | | Traceback correctness | 3 | ✓ | | Local alignment subset | 1 | ✓ | | Multiple alignments | 2 | ✓ | | Single character sequences | 2 | ✓ | | Reversibility (score symmetric) | 1 | ✓ | | Affine vs linear gap equivalence | 1 | ✓ | | Real protein example | 3 | ✓ | **Edge cases tested:** - **Empty sequences** (seq1 empty, seq2 empty, both empty) → score 0 - **No match possible** (AAAA vs CCCC with very high gap penalties) → score 0 - **Single characters** (A vs A = match, A vs B = mismatch) - **Affine gap penalties** (opening vs extension costs behave correctly) - **Reversibility** (swapping sequences yields identical score) - **Non-standard AA** (`'O'` raises `KeyError` as expected — not in BLOSUM62) - **`X` (any amino acid)** works correctly (X vs X = −1 in BLOSUM62) **Example output** from a canonical globin alignment: ``` Score: 90.0 Aln1: LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF Aln2: LTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRFFESF ``` ---
{"model": "pi", "problem_class": "smith-waterman-affine", "result": "passed", "tests": 24}