3: Local Sequence Alignment: Smith-Waterman Algorithm

Chapter 3: Local Sequence Alignment: Smith-Waterman Algorithm

While global alignment forces entire sequences to align from end to end, biology frequently presents conserved domains or active site motifs embedded within highly divergent non-homologous flanking sequences. Published in 1981 by Temple Smith and Michael Waterman, the Smith-Waterman algorithm identifies the optimal local region of similarity between two sequences.

1. Key Differences from Global Alignment

• Initialization: The first row and column are initialized to 0: \(H(i, 0) = 0\) and \(H(0, j) = 0\). Leading terminal gaps are never penalized.
• Traceback Initiation: Traceback does not start at the bottom-right corner. Instead, the algorithm scans the entire matrix to locate the cell containing the maximum value \(\max_{i,j} H(i, j)\).
• Traceback Termination: Traceback terminates as soon as a cell with a value of 0 is encountered.

2. Python Implementation

def smith_waterman(seq1, seq2, match=3, mismatch=-3, gap=-2):
    m, n = len(seq1), len(seq2)
    H = [[0] * (n + 1) for _ in range(m + 1)]
    max_score = 0
    max_pos = (0, 0)
    
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            score_diag = H[i-1][j-1] + (match if seq1[i-1] == seq2[j-1] else mismatch)
            score_up   = H[i-1][j] + gap
            score_left = H[i][j-1] + gap
            H[i][j] = max(0, score_diag, score_up, score_left)
            if H[i][j] > max_score:
                max_score = H[i][j]
                max_pos = (i, j)
                
    # Traceback from maximum position until 0
    align1, align2 = [], []
    i, j = max_pos
    while H[i][j] > 0:
        if H[i][j] == H[i-1][j-1] + (match if seq1[i-1] == seq2[j-1] else mismatch):
            align1.append(seq1[i-1])
            align2.append(seq2[j-1])
            i -= 1; j -= 1
        elif H[i][j] == H[i-1][j] + gap:
            align1.append(seq1[i-1])
            align2.append('-')
            i -= 1
        else:
            align1.append('-')
            align2.append(seq2[j-1])
            j -= 1
            
    return max_score, ''.join(reversed(align1)), ''.join(reversed(align2))

score, a1, a2 = smith_waterman("TGTTACGG", "GGTTGACTA")
print(f"Max Local Score: {score}\nLocal 1: {a1}\nLocal 2: {a2}")