dynamic_programming.needleman_wunsch

Needleman-Wunsch algorithm for global sequence alignment.

Reference:

https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm

The Needleman-Wunsch algorithm (1970) is a dynamic programming algorithm used in bioinformatics and computational biology to find the optimal global alignment between two sequences (such as DNA, RNA, or protein sequences).

Unlike local alignment algorithms (e.g., Smith-Waterman), which find the highest-scoring local sub-regions, Needleman-Wunsch aligns both sequences across their entire lengths from start to finish.

Algorithm: 1. Initialization:

  • Construct a matrix of size (m + 1) x (n + 1) where m and n are sequence lengths.

  • Initialize boundary conditions: score_matrix[i][0] = i * gap_score score_matrix[0][j] = j * gap_score

  1. Matrix Filling (Recurrence Relation): For each cell (i, j):

    diagonal = score_matrix[i - 1][j - 1] + (match_score if seq1[i-1] == seq2[j-1]

    else mismatch_score)

    deletion = score_matrix[i - 1][j] + gap_score insertion = score_matrix[i][j - 1] + gap_score score_matrix[i][j] = max(diagonal, deletion, insertion)

  2. Traceback: - Start from the bottom-right cell (m, n) and trace back to (0, 0). - At each step, determine which direction (diagonal, up, or left) produced the

    maximum score, assembling the aligned sequences in reverse order.

Complexity:

Time Complexity: O(m * n) where m and n are the lengths of the sequences. Space Complexity: O(m * n) to store the score matrix for traceback.

Functions

needleman_wunsch(→ tuple[str, str, int])

Compute the optimal global sequence alignment using Needleman-Wunsch.

Module Contents

dynamic_programming.needleman_wunsch.needleman_wunsch(sequence1: str, sequence2: str, match_score: int = 1, mismatch_score: int = -1, gap_score: int = -1) → tuple[str, str, int]

Compute the optimal global sequence alignment using Needleman-Wunsch.

Parameters:

sequence1: The first input sequence to align. sequence2: The second input sequence to align. match_score: Score awarded when two characters match (default: 1). mismatch_score: Penalty score when characters do not match (default: -1). gap_score: Penalty score for introducing a gap ‘-’ (default: -1).

Returns:

A tuple containing: - aligned_sequence1: The first aligned sequence with inserted gaps. - aligned_sequence2: The second aligned sequence with inserted gaps. - alignment_score: The total optimal alignment score.

Raises:

ValueError: If gap_score is positive (gap must be neutral or a penalty).

Examples:
>>> # Wikipedia classic example
>>> needleman_wunsch(
...     "GCATGCG", "GATTACA", match_score=1, mismatch_score=-1, gap_score=-1
... )
('GCA-TGCG', 'G-ATTACA', 0)

>>> # Identical sequences
>>> needleman_wunsch(
...     "ACGT", "ACGT", match_score=2, mismatch_score=-1, gap_score=-2
... )
('ACGT', 'ACGT', 8)

>>> # Completely mismatched sequences
>>> needleman_wunsch(
...     "AAAA", "TTTT", match_score=1, mismatch_score=-1, gap_score=-2
... )
('AAAA', 'TTTT', -4)

>>> # One sequence is empty
>>> needleman_wunsch("AGTC", "", match_score=1, mismatch_score=-1, gap_score=-1)
('AGTC', '----', -4)

>>> # Both sequences are empty
>>> needleman_wunsch("", "")
('', '', 0)

>>> # Protein sequence example
>>> needleman_wunsch(
...     "HEAGAWGHEE", "PAWHEAE", match_score=2, mismatch_score=-1, gap_score=-2
... )
('HEAGAWGHE-E', '---PAW-HEAE', -1)

>>> # Invalid gap score
>>> needleman_wunsch("A", "C", gap_score=5)
Traceback (most recent call last):
    ...
ValueError: gap_score must be non-positive (<= 0)