| 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"
hbb = "MVHLTPEEKSAVTALWGKVNVDE"
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)
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)
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)
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-objectsbio-alignment-msa-parsingbio-database-access-blast-searchesbio-alignment-alignment-io
