# Bio Core Pairwise Sequence Alignment

> Align two protein/DNA sequences with Biopython's PairwiseAligner (global Needleman-Wunsch, local Smith-Waterman), BLOSUM/PAM substitution matrices, and affine gap penalties; compute percent identity/similarity and E-values. Use when doing pairwise sequence alignment, ortholog comparison, dot plots, or choosing BLOSUM62 vs PAM250.

- Skill: `pavel-kravchenko/bio-core-pairwise-sequence-alignment` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-pairwise-sequence-alignment`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-pairwise-sequence-alignment/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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).

```python
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.

```python
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.

```python
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`