bio-core-pairwise-sequence-alignment

name	bio-core-pairwise-sequence-alignment
description	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.
tool_type	python
primary_tool	biopython

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 - ]  mismatch
            diag = F[i - , j - ] + s
            up = F[i - , j] + gap
            left = F[i, j - ] + gap
            F[i, j] = (diag, up, left)
            T[i, j] =   F[i, j] == diag  (  F[i, j] == up  )

    aligned1, aligned2 = [], []
    i, j = n, m
     i >   j > :
         i >   j >   T[i, j] == :
            aligned1.append(seq1[j - ]); aligned2.append(seq2[i - ]); i -= ; j -= 
         i >   T[i, j] == :
            aligned1.append(); aligned2.append(seq2[i - ]); i -= 
        :
            aligned1.append(seq1[j - ]); aligned2.append(); j -= 

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

a1, a2, score, F = needleman_wunsch(, , =, mismatch=-, gap=-)
(a1); (a2); (, 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] = 
             F[i, j] == up:
                T[i, j] = 
            :
                T[i, j] = 
             F[i, j] > max_score:
                max_score, max_i, max_j = F[i, j], i, j

    aligned1, aligned2 = [], []
    i, j = max_i, max_j
     i >   j >   T[i, j] != :
         T[i, j] == :
            aligned1.append(seq1[j - ]); aligned2.append(seq2[i - ]); i -= ; j -= 
         T[i, j] == :
            aligned1.append(); aligned2.append(seq2[i - ]); i -= 
        :
            aligned1.append(seq1[j - ]); aligned2.append(); j -= 

     .join((aligned1)), .join((aligned2)), max_score, F


 () -> :
    
     sub_matrix  :
         Bio.Align  substitution_matrices
        sub_matrix = substitution_matrices.load()
    identities = similarities = gaps = 
     a, b  (aligned1, aligned2):
         a ==   b == :
            gaps += 
         a == b:
            identities += ; similarities += 
         sub_matrix[a, b] > :
            similarities += 
    aligned_pos = (aligned1) - gaps
     {
        :  * identities / aligned_pos  aligned_pos  ,
        :  * similarities / aligned_pos  aligned_pos  ,
        :  * gaps / (aligned1)  aligned1  ,
    }

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