smith waterman algorithm choose more than one alignment

Beyond the Best: Finding Multiple Optimal Alignments with Smith-Waterman

Explore how to modify the Smith-Waterman algorithm to identify not just one, but all optimal local sequence alignments, crucial for comprehensive bioinformatics analysis.

The Smith-Waterman algorithm is a cornerstone in bioinformatics for identifying local similarities between two sequences. Unlike its global counterpart, Needleman-Wunsch, Smith-Waterman excels at finding conserved regions within otherwise divergent sequences. However, the standard implementation typically yields only one optimal alignment, even when multiple equally good alignments exist. This article delves into how to extend the Smith-Waterman algorithm to discover all such optimal local alignments, providing a more complete picture of sequence homology.

Understanding the Standard Smith-Waterman Algorithm

Before we can modify the algorithm to find multiple alignments, it's essential to grasp its core mechanics. Smith-Waterman uses dynamic programming to fill a scoring matrix, where each cell H(i, j) represents the optimal local alignment score ending at position i in sequence A and j in sequence B. The recurrence relation for H(i, j) is defined as:

H(i, j) = max(0, H(i-1, j-1) + S(A_i, B_j), H(i-1, j) - d, H(i, j-1) - d)

Where S(A_i, B_j) is the similarity score for aligning characters A_i and B_j, and d is the gap penalty. The max(0, ...) term ensures that local alignments start and end anywhere, and negative scores are reset to zero, effectively breaking off poor alignments. The highest value in the entire matrix indicates the score of the best local alignment.

flowchart TD
    A[Initialize Matrix H and P] --> B{Fill Matrix H}
    B --> C{Calculate H(i,j) using recurrence relation}
    C --> D{Track maximum score in H}
    D --> E{Find cell with maximum score}
    E --> F{Traceback from max score until H(i,j) = 0}
    F --> G[Output Single Optimal Alignment]

Standard Smith-Waterman Algorithm Flow

Modifying for Multiple Optimal Alignments

To find all optimal alignments, we need to adapt the traceback phase. The standard approach performs a single traceback from the highest scoring cell until a cell with a score of zero is reached. For multiple alignments, we must identify all cells that achieve the maximum score and initiate a traceback from each of them. Furthermore, during traceback, if a cell H(i, j) could have been reached from multiple preceding cells (e.g., H(i-1, j-1), H(i-1, j), or H(i, j-1)) with the same optimal score, we must explore all these paths. This implies that instead of storing a single pointer for traceback, we need to store a list of all possible preceding cells that led to the current optimal score.

Implementation Strategy: Python 2.7 Example

Let's outline a Python 2.7 implementation strategy. We'll need functions for matrix initialization, filling, and a recursive traceback function that can handle multiple paths. The traceback function will need to keep track of the current alignment and add it to a list of results once a zero score is encountered.

def smith_waterman_multiple(seq1, seq2, match=2, mismatch=-1, gap=-1):
    n, m = len(seq1), len(seq2)
    H = [[0] * (m + 1) for _ in range(n + 1)]
    P = [[[] for _ in range(m + 1)] for _ in range(n + 1)] # Pointer matrix for multiple paths
    max_score = 0
    max_pos = []

    for i in range(1, n + 1):
        for j in range(1, m + 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

            scores = [0, score_diag, score_up, score_left]
            H[i][j] = max(scores)

            # Store all directions that lead to the max score
            if H[i][j] == score_diag and H[i][j] != 0: P[i][j].append('D')
            if H[i][j] == score_up and H[i][j] != 0: P[i][j].append('U')
            if H[i][j] == score_left and H[i][j] != 0: P[i][j].append('L')

            if H[i][j] > max_score:
                max_score = H[i][j]
                max_pos = [(i, j)]
            elif H[i][j] == max_score and H[i][j] != 0:
                max_pos.append((i, j))

    alignments = []
    for r, c in max_pos:
        traceback_multiple(H, P, seq1, seq2, r, c, '', '', alignments)

    return max_score, alignments

def traceback_multiple(H, P, seq1, seq2, i, j, current_align1, current_align2, alignments):
    if H[i][j] == 0:
        alignments.append((current_align1[::-1], current_align2[::-1])) # Reverse to get correct order
        return

    for direction in P[i][j]:
        if direction == 'D':
            traceback_multiple(H, P, seq1, seq2, i-1, j-1, current_align1 + seq1[i-1], current_align2 + seq2[j-1], alignments)
        elif direction == 'U':
            traceback_multiple(H, P, seq1, seq2, i-1, j, current_align1 + seq1[i-1], current_align2 + '-', alignments)
        elif direction == 'L':
            traceback_multiple(H, P, seq1, seq2, i, j-1, current_align1 + '-', current_align2 + seq2[j-1], alignments)

# Example Usage:
# seq1 = "GGATCGA"
# seq2 = "GAATTC"
# score, alns = smith_waterman_multiple(seq1, seq2)
# print "Max Score:", score
# for a1, a2 in alns:
#     print a1
#     print a2
#     print

Python 2.7 code for Smith-Waterman with multiple alignment traceback.

Analyzing the Results and Further Considerations

The smith_waterman_multiple function will return the maximum score and a list of all unique optimal alignments. Each alignment in the list will be a tuple of two strings representing the aligned segments of seq1 and seq2. It's important to note that the definition of 'optimal' here refers to alignments achieving the highest possible score. If you need to find suboptimal alignments (e.g., alignments with scores within a certain percentage of the maximum), the algorithm would require further modifications, potentially involving a more complex traceback or iterative removal of already found optimal paths.

This extended Smith-Waterman algorithm is particularly useful in scenarios where subtle variations in alignment can have significant biological implications, such as identifying alternative splicing sites, protein domain boundaries, or conserved regulatory elements that might be missed by a single-best alignment approach.

flowchart TD
    A[Initialize Matrix H and P] --> B{Fill Matrix H}
    B --> C{Calculate H(i,j) and store ALL optimal directions in P[i][j]}
    C --> D{Track ALL cells with maximum score}
    D --> E{For EACH max score cell (r,c)}
    E --> F{Initiate recursive traceback_multiple(H, P, ..., r, c)}
    F --> G{Explore ALL paths from (r,c) until H(i,j) = 0}
    G --> H[Collect ALL unique optimal alignments]

Modified Smith-Waterman Algorithm Flow for Multiple Alignments