Rank Alignment – Accurately Overlap Two DNAs Using Dynamic Programming

After Completing This Topic

By combining the dynamic programming, recursion-vs-iteration, and 2D matrix concepts learned in the textbook, you will be able to create your own alignment tool that accurately aligns two DNA sequences and identifies the locations of insertions, deletions, and substitutions. You will gain a deeper understanding of the underlying principles of this alignment process, which BLAST performs seemingly magically.

This article is an educational general example. Mutation analysis is a common task in bioinformatics, so it was chosen as the subject matter.

"Why Don't the Mutation Positions Match?" — The Pitfall of Naive Comparison

Let's say you're comparing a wild-type sequence of the BRCA1 gene with a mutated sequence obtained from a patient.

Wild-type:   ATGCTAGCATGCA
Mutated: ATGCAGCATGCA

The simplest comparison is to match them one position at a time.

def compare_naive(seq1: str, seq2: str) -> list[tuple[int, str, str]]:
    diffs = []
    for i in range(min(len(seq1), len(seq2))):
        if seq1[i] != seq2[i]:
            diffs.append((i, seq1[i], seq2[i]))
    return diffs

If you run this on the two sequences above, you'll get something like this:

Position 4: T → A
Position 5: A → G
Position 6: G → C
Position 7: C → A
Position 8: A → T
Position 9: T → G
Position 10: G → C
Position 11: C → A

Strange, isn't it? After position 4, almost every position is different. Is there really that many mutations?

In reality, if you look at the two sequences, the answer is obvious. In the wild-type sequence, the T at position 4 has been deleted.

Wild-type:   ATGC[T]AGCATGCA
Mutated: ATGC AGCATGCA
                ↑ The offset shifts from here

A single letter deletion causes all subsequent positions to shift, making it look like there are multiple substitutions. A naive comparison doesn't notice this deletion or insertion.

The solution to this problem is sequence alignment. It involves inserting as many gaps (represented by the - character) as needed between the two sequences to correct for the misalignment.

Wild-type:   ATGCTAGCATGCA
Mutated: ATGC-AGCATGCA   ← Gap at position 4

Now, after position 5, everything matches perfectly. It becomes clear that the true mutation is a single deletion.

The problem is deciding where to put these gaps. If the two sequences are about 100 characters long, the number of possible locations to insert gaps increases exponentially. How do you find the combination that represents the "most natural" alignment? This article builds the answer from the ground up.

Let's Start by Looking at the Finished Product (Run the Black Box First)

The tool we're going to build takes two sequences and returns the optimal alignment.

aligned1, aligned2, score = align(seq1, seq2)

print(aligned1)  # 'ATGCTAGCATGCA'
print(aligned2)  # 'ATGC-AGCATGCA'
print(score)     # 8 (example)

=== Alignment Result ===
Wild:   ATGCTAGCATGCA
Mut:    ATGC-AGCATGCA
Score:  8

=== Naive Comparison Without Alignment ===
Wild:   ATGCTAGCATGCA
Mut:    ATGCAGCATGCA
Mismatched Positions: 8 (almost all)

=== Actual Variation After Alignment ===
Position 4: T deletion — 1 true variation

These are the same two sequences, but the interpretation changes completely before and after alignment. Alignment is the first and crucial step in variation analysis. Now, let's explore the principles behind creating this alignment.

What components are used to assemble this tool (component breakdown)?

Sequence Aligner
   ┌──────────────────────────────────────────────────┐
   │  [Input] Load two sequences ─── Component: String/File I/O  │  ← Provided as a complete tool
   │              │                                      │
   │              ▼                                      │
   │  [Step 1] Fill the score matrix                          │
   │        Component: 2D Array + Dynamic Programming         │  ← To be implemented from scratch ★
   │              │                                      │
   │              ▼                                      │
   │  [Step 2] Restore the optimal path by backtracking       │
   │        Component: Recursion vs. Iteration               │  ← To be implemented from scratch ★
   │              │                                      │
   │              ▼                                      │
   │  [Output] Aligned sequences + score                      │
   └──────────────────────────────────────────────────┘

📌 If you are unfamiliar with these concepts (links at the top)

• Dynamic Programming — Avoiding exponential complexity by storing solutions to subproblems
• Recursion vs. Iteration — A trade-off between elegance and efficiency
• 2D Array — Values stored in a grid

The new concepts to be implemented from scratch are Dynamic Programming, 2D Matrix, and Recursion/Iteration Selection. Sequence handling is provided as a complete tool. There are only three of them, which is within cognitive limits.

Step 1: Preparing the Sequences (Provided)

First, we prepare two sequences to be aligned. In a real application, these would be read from a FASTA file, but for this browser-based exercise, we define them as strings.

# Wild-type BRCA1 partial sequence (for demonstration purposes)
wild = "ATGCTAGCATGCA"

# Mutant sequence (with a T deletion at position 4)
mut = "ATGCAGCATGCA"

assert len(wild) == 13
assert len(mut) == 12
assert all(c in "ACGT" for c in wild)
assert all(c in "ACGT" for c in mut)

Our goal is to align these two sequences by inserting gaps (-) appropriately.

Step 2: Defining Scoring Rules (Fully Implemented)

First, we need to define how to measure the "quality" of an alignment. A standard approach is to use a scoring system.

MATCH = 1        # +1 if the same bases are aligned
MISMATCH = -1    # -1 if different bases are aligned
GAP = -2         # -2 for each gap

def score_pair(a: str, b: str) -> int:
    """The alignment score for two bases (or a gap)."""
    if a == b:
        return MATCH
    return MISMATCH

assert score_pair("A", "A") == 1
assert score_pair("A", "T") == -1

This scoring rule determines the "preference" of the alignment algorithm.

• If the gap penalty is low (e.g., -0.5), the algorithm will prefer alignments with many gaps.
• If the gap penalty is high (e.g., -5), the algorithm will prefer mismatches over gaps.

In practice, more sophisticated scoring matrices like BLOSUM/PAM are used. Here, we use a simple rule for learning purposes.

Step 3: Build It — Start with Recursion, Then Experience the Explosion ★ (Recursion vs. Iteration)

✍️ Fill-in section. Component = Recursion. First, implement the alignment score calculation using recursion, and then experience why it's not practical.

Expressing the alignment score problem recursively is surprisingly simple. Starting from the ends of the two sequences, we recursively try three choices at each position.

• Align the last characters of the two sequences → score_pair(s1[-1], s2[-1]) + align(s1[:-1], s2[:-1])
• Align the last character of s1 with a gap → GAP + align(s1[:-1], s2)
• Align the last character of s2 with a gap → GAP + align(s1, s2[:-1])

We choose the option with the highest score out of these three. Writing this recursively gives us:

def align_recursive(s1: str, s2: str) -> int:
    """Recursively calculate only the score of the optimal alignment (no traceback)."""
    if not s1:
        return len(s2) * GAP
    if not s2:
        return len(s1) * GAP
    return max(
        align_recursive(s1[:-1], s2[:-1]) + score_pair(s1[-1], s2[-1]),
        align_recursive(s1[:-1], s2) + GAP,
        align_recursive(s1, s2[:-1]) + GAP,
    )

# Verification: Only with short examples
assert align_recursive("AA", "AA") == 2       # Two matches
assert align_recursive("AT", "A") == 1 + GAP  # Match + gap = 1 - 2 = -1

This does give the correct answer. However, it's completely unusable for long sequences. If you run the following code, you'll see why.

import time

# Even with a length of around 15, it becomes a recursive nightmare
short1 = "ATGCTAGCATGCAAG"
short2 = "ATGCAGCATGCAAG"

t0 = time.time()
score = align_recursive(short1, short2)
elapsed = time.time() - t0
print(f"Length 15 Recursion: {elapsed:.2f} seconds")
# It takes several seconds. For length 20, it takes minutes, and for 25, it takes hours.

The reason is that recursion repeatedly calculates the same subproblems hundreds of times. The value of align_recursive("ATGCT", "ATGC") is calculated repeatedly in multiple branches of the recursion tree, each time from the beginning.

This is the power of O(3^n) time complexity. At this point, dynamic programming comes into play.

🤔 Self-explanatory prompt Draw the recursion tree for align_recursive("AT", "A"). How many subproblems appear repeatedly? (Even expanding to lengths of 3 and 4 makes the overlaps obvious.)

Step 4: Storing Answers in a Grid ★ (2D Matrix + Dynamic Programming)

✍️ The core concept: fill in the grid. The components are: 2D array + dynamic programming. Store answers to subproblems in the grid to avoid redundant calculations.

The idea is this: store the answer to align(s1[:i], s2[:j]) in the grid dp[i][j]. This way, we don't need to compute the same subproblem twice.

🔎 What is Dynamic Programming? (Drawer — dynamic-programming) When the solution to a large problem can be constructed from the solutions to smaller problems, store the solutions to the smaller problems so you don't have to compute them twice. This turns the exponential explosion of recursion into polynomial time with a single memoization. Filling in a grid is specifically called "Bottom-Up Dynamic Programming."

The grid size is (len(s1)+1) × (len(s2)+1). The extra row and column are for the base case: aligning with the empty string.

def build_score_matrix(s1: str, s2: str) -> list[list[int]]:
    """dp[i][j] = optimal alignment score for s1[:i] and s2[:j]."""
    n, m = len(s1), len(s2)
    dp = [[0] * (m + 1) for _ in range(n + 1)]

    # Base case: aligning with an empty string involves only gaps
    for i in range(1, n + 1):
        dp[i][0] = dp[i-1][0] + GAP
    for j in range(1, m + 1):
        dp[0][j] = dp[0][j-1] + GAP

    # Fill in the grid (Bottom-Up)
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            match = dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1])
            delete = dp[i-1][j] + GAP  # s1 character and a gap
            insert = dp[i][j-1] + GAP  # gap and s2 character
            dp[i][j] = max(match, delete, insert)
    return dp

dp = build_score_matrix(wild, mut)

# Verification: the final cell should match the recursive answer
assert dp[len(wild)][len(mut)] == align_recursive(wild, mut)
# Size verification
assert len(dp) == len(wild) + 1
assert len(dp[0]) == len(mut) + 1
# Leftmost column should be a pure accumulation of gaps
assert dp[3][0] == 3 * GAP

Get the same answer, but much faster.

import time

# A string of length 15, which took a few seconds with recursion
t0 = time.time()
score = build_score_matrix(short1, short2)[len(short1)][len(short2)]
elapsed = time.time() - t0
print(f"Length 15 dynamic programming: {elapsed*1000:.2f}ms")
# In milliseconds

What took a few seconds with recursion now takes milliseconds. This is the power of O(n×m). Each cell in the grid is computed only once.

🔎 Why is a Grid the Answer? (Drawer — matrix-2d-array) A subproblem that depends on two indices (i, j) is naturally represented as a 2D grid. Each cell is a subproblem, and the arrows between cells represent the dependencies. If you can visualize this structure, it becomes easier to approach DP problems.

🤔 Self-Explanatory Prompt When computing dp[i][j], you only need dp[i-1][j-1], dp[i-1][j], and dp[i][j-1]. Do you really need to store the entire grid? In fact, you only need to store the previous row to fill in the next row, so you can compute the score with O(min(n,m)) space. When would this optimization not be possible? (Hint: traceback)

Step 5: Reconstructing the Alignment by Tracing Back the Grid ★

✍️ The Fill-In-The-Blanks Section. The component = Iterative Backtracking. Start from the last cell of the grid and trace back how the value was created.

Knowing only the score is different from knowing exactly how the alignment was made. To reconstruct the actual alignment (where the gaps are inserted), we need to trace back the grid from the end to the beginning.

In each cell dp[i][j], we trace back which of the three candidates we chose.

• If we chose dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]) → move diagonally (match/mismatch)
• If we chose dp[i-1][j] + GAP → move up (s1 character and gap)
• If we chose dp[i][j-1] + GAP → move left (gap and s2 character)

def traceback(dp: list[list[int]], s1: str, s2: str) -> tuple[str, str]:
    aligned1: list[str] = []
    aligned2: list[str] = []
    i, j = len(s1), len(s2)

    while i > 0 or j > 0:
        current = dp[i][j]
        if i > 0 and j > 0 and current == dp[i-1][j-1] + score_pair(s1[i-1], s2[j-1]):
            aligned1.append(s1[i-1])
            aligned2.append(s2[j-1])
            i -= 1
            j -= 1
        elif i > 0 and current == dp[i-1][j] + GAP:
            aligned1.append(s1[i-1])
            aligned2.append("-")
            i -= 1
        else:
            aligned1.append("-")
            aligned2.append(s2[j-1])
            j -= 1

    return "".join(reversed(aligned1)), "".join(reversed(aligned2))

a1, a2 = traceback(dp, wild, mut)

# Verification: Length is the same, original sequences are restored when gaps are removed
assert len(a1) == len(a2)
assert a1.replace("-", "") == wild
assert a2.replace("-", "") == mut
# The gap we inserted should actually be found
assert "-" in a2  # There should be a gap in the mutant sequence

This function is implemented with iteration, but the same logic can also be written recursively. When the grid is very large, iteration is better for preventing stack overflow, but recursion can sometimes be more readable. The choice depends on the problem size.

🔎 When to Use Recursion or Iteration (Drawer — recursion-vs-iteration) Recursion is elegant when it naturally expresses the problem. Iteration is safer because it doesn't stack up and is often faster. Small problem, elegance needed → recursion / Large problem, stability needed → iteration. In this tool, step 3 (filling the grid) is implemented with iteration, and step 4 (traceback) is also implemented with iteration. Both can be implemented with recursion, but we choose iteration because of the risk of stack overflow when the sequence becomes long.

🤔 Self-Explanatory Prompt In traceback, we choose one of three candidates. If the scores of two candidates are the same (tie), which one should we choose? Could it affect the actual alignment result? (Hint: There may be multiple optimal alignments.)

Combining the Pieces: The Complete Aligner Class

Let's bundle the two functions into a single tool.

class Aligner:
    def __init__(self, match=1, mismatch=-1, gap=-2):
        self.match = match
        self.mismatch = mismatch
        self.gap = gap

    def _score(self, a: str, b: str) -> int:
        return self.match if a == b else self.mismatch

    def align(self, s1: str, s2: str) -> tuple[str, str, int]:
        n, m = len(s1), len(s2)
        dp = [[0] * (m + 1) for _ in range(n + 1)]
        for i in range(1, n + 1):
            dp[i][0] = dp[i-1][0] + self.gap
        for j in range(1, m + 1):
            dp[0][j] = dp[0][j-1] + self.gap
        for i in range(1, n + 1):
            for j in range(1, m + 1):
                dp[i][j] = max(
                    dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]),
                    dp[i-1][j] + self.gap,
                    dp[i][j-1] + self.gap,
                )
        # Backtracking
        a1: list[str] = []
        a2: list[str] = []
        i, j = n, m
        while i > 0 or j > 0:
            if i > 0 and j > 0 and dp[i][j] == dp[i-1][j-1] + self._score(s1[i-1], s2[j-1]):
                a1.append(s1[i-1]); a2.append(s2[j-1]); i -= 1; j -= 1
            elif i > 0 and dp[i][j] == dp[i-1][j] + self.gap:
                a1.append(s1[i-1]); a2.append("-"); i -= 1
            else:
                a1.append("-"); a2.append(s2[j-1]); j -= 1
        return "".join(reversed(a1)), "".join(reversed(a2)), dp[n][m]


aligner = Aligner()
a1, a2, score = aligner.align(wild, mut)
print(a1)
print(a2)
print(f"score = {score}")

# Verification: the class result matches the function result
dp = build_score_matrix(wild, mut)
assert score == dp[len(wild)][len(mut)]

This is a miniature version of the Needleman-Wunsch algorithm, which is what we've been discussing today. Real-world tools (like BLAST or needle from EMBOSS) simply add more sophisticated scoring matrices, separate scores for gap opening and extension, local alignment options, and so on. But the underlying framework is the same as what you just created.

Performance Deep Dive — Why is it so fast?

import time

def bench(n: int, m: int):
    s1 = "AT" * (n // 2)
    s2 = "AG" * (m // 2)
    t0 = time.time()
    build_score_matrix(s1, s2)
    return time.time() - t0

for size in [10, 50, 100, 200]:
    t = bench(size, size)
    print(f"Length {size:4d}: {t*1000:.2f}ms  ({size*size} cells)")
assert bench(100, 100) < bench(200, 200)

The numbers vary from machine to machine, but the trend is always the same. It is directly proportional to the number of cells (n×m). The critical difference is that while the recursive version explodes at 3^n, the DP version grows at n×m.

🔎 Summarized in Big-O (Drawer — big-o-notation)

• Naive recursive version: O(3^n) — Exponential. Impractical for lengths of 30.
• Dynamic programming version: O(n·m) — Polynomial. Can handle lengths of 10,000 in a few seconds.
• Storing (in a grid) one thing turned the exponential into a polynomial. This is the magic of dynamic programming.

There Are Other Paths (Multipass Reflection)

• Global vs. Local Alignment: What we created is Needleman-Wunsch, which aligns the entire sequence from beginning to end (global). In practice, Smith-Waterman, which finds fragments of one sequence within another, is often used (local). The algorithmic framework is almost the same, but only the scoring system and initial conditions differ. When to use which: For aligning evolutionary sequences = global / For finding domains within a gene = local.
• BLAST's Heuristic Approach: Our dynamic programming's complexity is proportional to the lengths of the two sequences, n·m, which is still large for genomes (billions of base pairs). BLAST sacrifices exact solutions in favor of quickly finding approximate solutions using a heuristic (seed → extend). Trade-off: 100% optimal vs. practical speed.
• Refined Scoring Matrix: We use a binary match/mismatch rule, but in practice, we use matrices like BLOSUM62/PAM250, which assign different scores to each residue pair. This reflects the chemical similarity of protein sequences.
• Affine Gap Penalty: We use a fixed penalty of -2 for each gap, but in reality, we use a dual rule where the cost of opening a gap is high, but the cost of extending it is low (e.g., opening -10, extension -1). This reflects the evolutionary observation that a single long gap is more natural than multiple short gaps.

Key takeaway: "Store possible subproblems in a grid to avoid redundant calculations." This principle can be reused in countless applications beyond sequence alignment, such as pathfinding, edit distance, and knowledge retrieval. The framework you just created opens up many possibilities.

Next Steps (Links to Further Exploration)

• Why do we visit the cells in this order when filling the grid? → Dependencies and Visiting Order
• Memory optimization when the sequence is very long → 1D DP Compression Technique
• Combining it with the index we created earlier → Application: Sequence Database Indexing

Try It Yourself (Independent Exercises)

1. Insertion Scenario: Align the wild-type sequence ATGCATGCA with the mutant sequence ATGCXATGCA (where X represents an inserted base). Assert that the insertion position is correctly identified.
2. Score Tuning: How do the results change when using Aligner(gap=-5)? What about gap=-0.5? Summarize how the gap penalty determines the "character" of the alignment.
3. Multiple Optimal Alignments: Create a pair of sequences that have two optimal alignments. Extend the traceback function to create a version that returns all optimal alignments.
4. Challenge — Memory O(min(n,m)): If you only need the score, you can keep only the previous row instead of the entire grid. Implement this optimized version. (You will have to give up traceback, but it will be useful for very long sequences.)

Summary

We have conquered the problem of "finding the optimal alignment of two sequences" by breaking it down into three components:

• Dynamic programming transformed the exponential explosion (O(3^n)) into polynomial time (O(n·m)).
• A 2D matrix became the grid for storing the solutions to subproblems.
• The choice between recursion and iteration ensured stability in both filling the grid and backtracking.

When BLAST seemingly magically throws out alignment results, you can now see what's going on inside. It's not magic, but rather simply applying sequences to the dp grid learned in textbooks.

This article is a general educational example. Practical sequence alignment tools (BLAST/EMBOSS/Bowtie, etc.) add sophisticated scoring matrices, gap double penalties, seed-and-extend heuristics, and parallelization to this. You can either build upon this framework to create a detailed version or rely on validated tools.