2: Global Sequence Alignment: Needleman-Wunsch Algorithm

Chapter 2: Global Sequence Alignment: Needleman-Wunsch Algorithm

Sequence alignment is the most fundamental algorithmic procedure in computational biology. Published in 1970 by Saul Needleman and Christian Wunsch, the Needleman-Wunsch algorithm applied Richard Bellman's principle of dynamic programming to solve the global sequence alignment problem with guaranteed mathematical optimality.

1. Matrix Initialization and Fill

A global alignment forces both sequences to align end-to-end from their first to last characters:

• Initialization: The top-left corner is set to \(F(0, 0) = 0\). The first column is initialized as \(F(i, 0) = -i imes d\), and the first row as \(F(0, j) = -j imes d\). This penalizes terminal end gaps.
• Fill: The matrix is populated row-by-row or column-by-column, evaluating the maximum of the diagonal (match/mismatch), vertical (deletion / gap in \(S_2\)), and horizontal (insertion / gap in \(S_1\)) transitions.
• Traceback: Initiates at the bottom-right cell \(F(m, n)\) and follows pointer arrows back to \(F(0, 0)\) to reconstruct the optimal alignment path.

2. Pure Python Implementation

def needleman_wunsch(seq1, seq2, match=2, mismatch=-1, gap=-2):
    m, n = len(seq1), len(seq2)
    # Initialize DP score matrix and traceback matrix
    dp = [[0] * (n + 1) for _ in range(m + 1)]
    
    for i in range(m + 1):
        dp[i][0] = i * gap
    for j in range(n + 1):
        dp[0][j] = j * gap
        
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            score_diag = dp[i-1][j-1] + (match if seq1[i-1] == seq2[j-1] else mismatch)
            score_up   = dp[i-1][j] + gap
            score_left = dp[i][j-1] + gap
            dp[i][j] = max(score_diag, score_up, score_left)
            
    # Traceback from bottom-right to top-left
    align1, align2 = [], []
    i, j = m, n
    while i > 0 or j > 0:
        if i > 0 and j > 0 and dp[i][j] == dp[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 i > 0 and dp[i][j] == dp[i-1][j] + gap:
            align1.append(seq1[i-1])
            align2.append('-')
            i -= 1
        else:
            align1.append('-')
            align2.append(seq2[j-1])
            j -= 1
            
    return dp[m][n], ''.join(reversed(align1)), ''.join(reversed(align2))

# Test alignment on two DNA fragments
score, a1, a2 = needleman_wunsch("GATTACA", "GCATGCU")
print(f"Optimal Score: {score}\nAlign 1: {a1}\nAlign 2: {a2}")