Bio Core Pairwise Sequence Alignment

Pairwise Sequence Alignment

When to Use

• Aligning two full-length orthologous proteins or genes end-to-end (global alignment).
• Finding a shared domain or motif inside a longer sequence (local alignment).
• Choosing a substitution matrix (BLOSUM vs PAM) or gap penalty scheme (linear vs affine) for an alignment.
• Explaining/implementing Needleman-Wunsch or Smith-Waterman dynamic programming from scratch (teaching, debugging aligner output).
• Computing percent identity/similarity, or judging whether a hit is statistically significant (E-value, twilight zone).

Version Compatibility

Biopython >= 1.80 (Bio.Align.PairwiseAligner, Bio.Align.substitution_matrices), Python >= 3.9, NumPy >= 1.24, matplotlib >= 3.7.

Prerequisites

• pip install biopython numpy matplotlib
• Familiarity with Seq objects (see bio-sequence-manipulation-seq-objects).

Key Concepts

• Global vs. local: Needleman-Wunsch aligns end-to-end — use for full-length orthologs of similar length. Smith-Waterman finds the best-matching subregion — use for shared domains or short query vs. long subject. BLAST uses heuristic local alignment.
• BLOSUM numbering: higher number = built from more similar sequences (BLOSUM80 for close relatives, BLOSUM45 for distant, BLOSUM62 is the BLAST/Biopython default). PAM is inverted: higher PAM = more divergence modeled (PAM250 ≈ 250 accepted mutations per 100 residues, PAM1 = 1% divergence extrapolated by matrix exponentiation).
• Affine gap penalties: linear cost is d * k (same cost per gap position); affine is d + (k-1) * e with e < d — expensive to open a gap, cheap to extend it. Biologically realistic since indels occur in runs (e.g. replication slippage), not as scattered single-position gaps. Typical values: linear d=8; affine d=10, e=0.5.
• Score comparability: a raw score isn't comparable across alignments of different length/composition — use percent identity, and for database searches use bit score/E-value (E-value scales with database size; bit score does not).

Goal: quick global/local alignment with Biopython

Approach: PairwiseAligner handles both modes; load a named substitution matrix and set affine gap scores (negative = penalty).

from Bio.Align import PairwiseAligner, substitution_matrices
from Bio.Seq import Seq

blosum62 = substitution_matrices.load("BLOSUM62")

def align_pair(seq1: str, seq2: str, mode: str = "global",
               matrix=blosum62, open_gap: float = -10, extend_gap: float = -0.5):
    """Align two sequences with Biopython's PairwiseAligner.

    mode: 'global' (Needleman-Wunsch) or 'local' (Smith-Waterman).
    Returns the best-scoring Alignment object; use str(alignment) to print it.
    """
    aligner = PairwiseAligner()
    aligner.mode = mode
    aligner.substitution_matrix = matrix
    aligner.open_gap_score = open_gap
    aligner.extend_gap_score = extend_gap
    alignments = aligner.align(seq1, seq2)
    return alignments[0]

hba = "MVLSPADKTNVKAAWGKVGAHAG"  # human hemoglobin alpha N-term
hbb = "MVHLTPEEKSAVTALWGKVNVDE"  # human hemoglobin beta N-term
best = align_pair(hba, hbb, mode="global")
print(best)
print(f"Score: {best.score:.1f}")

motif = "AWGKVGAHAG"
target = "MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSH"
hit = align_pair(target, motif, mode="local")
print(hit)

Goal: implement Needleman-Wunsch (global) from scratch

Approach: fill an (n+1) x (m+1) DP matrix with the recurrence F(i,j) = max(diag + s, up + gap, left + gap), then trace back from the bottom-right corner. O(n*m) time and space.

import numpy as np

def needleman_wunsch(seq1: str, seq2: str, match: int = 1, mismatch: int = -1, gap: int = -2):
    """Global alignment via Needleman-Wunsch with a linear gap penalty.

    Returns (aligned1, aligned2, score, F) where F is the DP score matrix.
    """
    m, n = len(seq1), len(seq2)
    F = np.zeros((n + 1, m + 1), dtype=int)
    T = np.zeros((n + 1, m + 1), dtype=int)  # 0=diag, 1=up, 2=left

    for i in range(1, n + 1):
        F[i, 0] = gap * i
        T[i, 0] = 1
    for j in range(1, m + 1):
        F[0, j] = gap * j
        T[0, j] = 2

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            s = match if seq1[j - 1] == seq2[i - 1] else mismatch
            diag = F[i - 1, j - 1] + s
            up = F[i - 1, j] + gap
            left = F[i, j - 1] + gap
            F[i, j] = max(diag, up, left)
            T[i, j] = 0 if F[i, j] == diag else (1 if F[i, j] == up else 2)

    aligned1, aligned2 = [], []
    i, j = n, m
    while i > 0 or j > 0:
        if i > 0 and j > 0 and T[i, j] == 0:
            aligned1.append(seq1[j - 1]); aligned2.append(seq2[i - 1]); i -= 1; j -= 1
        elif i > 0 and T[i, j] == 1:
            aligned1.append("-"); aligned2.append(seq2[i - 1]); i -= 1
        else:
            aligned1.append(seq1[j - 1]); aligned2.append("-"); j -= 1

    return "".join(reversed(aligned1)), "".join(reversed(aligned2)), F[n, m], F

a1, a2, score, F = needleman_wunsch("HEAGAWGHEE", "PAWHEAE", match=2, mismatch=-1, gap=-2)
print(a1); print(a2); print("Score:", score)

Goal: implement Smith-Waterman (local) and score an alignment

Approach: same recurrence as Needleman-Wunsch but floor every cell at 0 (max(0, diag, up, left)); traceback starts at the max-scoring cell and stops at the first 0, isolating the best local sub-region.

def smith_waterman(seq1: str, seq2: str, match: int = 2, mismatch: int = -1, gap: int = -1):
    """Local alignment via Smith-Waterman with a linear gap penalty."""
    m, n = len(seq1), len(seq2)
    F = np.zeros((n + 1, m + 1), dtype=int)
    T = np.zeros((n + 1, m + 1), dtype=int)  # 0=stop, 1=diag, 2=up, 3=left
    max_score, max_i, max_j = 0, 0, 0

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            s = match if seq1[j - 1] == seq2[i - 1] else mismatch
            diag = F[i - 1, j - 1] + s
            up = F[i - 1, j] + gap
            left = F[i, j - 1] + gap
            F[i, j] = max(0, diag, up, left)  # zero floor = key difference from NW
            if F[i, j] == 0:
                T[i, j] = 0
            elif F[i, j] == diag:
                T[i, j] = 1
            elif F[i, j] == up:
                T[i, j] = 2
            else:
                T[i, j] = 3
            if F[i, j] > max_score:
                max_score, max_i, max_j = F[i, j], i, j

    aligned1, aligned2 = [], []
    i, j = max_i, max_j
    while i > 0 and j > 0 and T[i, j] != 0:
        if T[i, j] == 1:
            aligned1.append(seq1[j - 1]); aligned2.append(seq2[i - 1]); i -= 1; j -= 1
        elif T[i, j] == 2:
            aligned1.append("-"); aligned2.append(seq2[i - 1]); i -= 1
        else:
            aligned1.append(seq1[j - 1]); aligned2.append("-"); j -= 1

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


def alignment_statistics(aligned1: str, aligned2: str, sub_matrix=None) -> dict:
    """Percent identity/similarity and gap fraction for an already-aligned pair."""
    if sub_matrix is None:
        from Bio.Align import substitution_matrices
        sub_matrix = substitution_matrices.load("BLOSUM62")
    identities = similarities = gaps = 0
    for a, b in zip(aligned1, aligned2):
        if a == "-" or b == "-":
            gaps += 1
        elif a == b:
            identities += 1; similarities += 1
        elif sub_matrix[a, b] > 0:
            similarities += 1
    aligned_pos = len(aligned1) - gaps
    return {
        "identity_pct": 100 * identities / aligned_pos if aligned_pos else 0,
        "similarity_pct": 100 * similarities / aligned_pos if aligned_pos else 0,
        "gap_pct": 100 * gaps / len(aligned1) if aligned1 else 0,
    }

E-value scales with database size (E ~= m * n * 2**(-bit_score)); the same alignment looks "more significant" in a smaller database, so always compare bit scores (database-independent) when judging hits across searches. Percent identity below ~20% is the "midnight zone" (indistinguishable from chance); 20-30% is the "twilight zone" — check E-value and alignment length/coverage before calling sequences homologous.

Pitfalls

• BLOSUM/PAM numbering is easy to get backwards — BLOSUM80 is for close sequences, PAM250 is for distant ones.
• A high raw score on a short alignment can outscore a biologically better long alignment — always report percent identity/coverage alongside the score.
• Global alignment (Needleman-Wunsch) on sequences of very different length forces spurious gaps across the whole unaligned region; use local (Smith-Waterman) or don't force full-length alignment when only a domain is shared.
• Pure-Python DP is O(n*m) in time and memory — fine for genes/proteins, but don't run it on whole chromosomes; use BLAST/minimap2 for genome-scale search instead.
• PairwiseAligner.align() can return many co-optimal alignments; alignments[0] is only one of possibly several equally-scoring paths.

See Also

• bio-sequence-manipulation-seq-objects
• bio-alignment-msa-parsing
• bio-database-access-blast-searches
• bio-alignment-alignment-io