Pairwise Sequence Alignment

Align two sequences using dynamic programming algorithms (Needleman-Wunsch for global, Smith-Waterman for local).

Required Import

from Bio.Align import PairwiseAligner from Bio.Seq import Seq from Bio import SeqIO

Core Concepts

Mode Algorithm Use Case

global

Needleman-Wunsch Full-length alignment, similar-length sequences

local

Smith-Waterman Find best matching regions, different-length sequences

Creating an Aligner

Basic aligner with defaults

aligner = PairwiseAligner()

Configure mode and scoring

aligner = PairwiseAligner(mode='global', match_score=2, mismatch_score=-1, open_gap_score=-10, extend_gap_score=-0.5)

For protein alignment with substitution matrix

from Bio.Align import substitution_matrices aligner = PairwiseAligner(mode='global', substitution_matrix=substitution_matrices.load('BLOSUM62'))

Performing Alignments

seq1 = Seq('ACCGGTAACGTAG') seq2 = Seq('ACCGTTAACGAAG')

Get all optimal alignments

alignments = aligner.align(seq1, seq2) print(f'Found {len(alignments)} optimal alignments') print(alignments[0]) # Print first alignment

Get score only (faster for large sequences)

score = aligner.score(seq1, seq2)

Alignment Output Format

target 0 ACCGGTAACGTAG 13 0 |||||.||||.|| 13 query 0 ACCGTTAACGAAG 13

Accessing Alignment Data

alignment = alignments[0]

Basic properties

print(alignment.score) # Alignment score print(alignment.shape) # (num_seqs, alignment_length) print(len(alignment)) # Alignment length

Get aligned sequences with gaps

target_aligned = alignment[0, :] # First sequence (target) with gaps query_aligned = alignment[1, :] # Second sequence (query) with gaps

Get coordinate mapping

print(alignment.aligned) # Array of aligned segment coordinates print(alignment.coordinates) # Full coordinate array

Alignment Counts (Identities, Mismatches, Gaps)

alignment = alignments[0] counts = alignment.counts()

print(f'Identities: {counts.identities}') print(f'Mismatches: {counts.mismatches}') print(f'Gaps: {counts.gaps}')

Calculate percent identity

total_aligned = counts.identities + counts.mismatches percent_identity = counts.identities / total_aligned * 100 print(f'Percent identity: {percent_identity:.1f}%')

Common Scoring Configurations

DNA/RNA Alignment

aligner = PairwiseAligner(mode='global', match_score=2, mismatch_score=-1, open_gap_score=-10, extend_gap_score=-0.5)

Protein Alignment

from Bio.Align import substitution_matrices blosum62 = substitution_matrices.load('BLOSUM62') aligner = PairwiseAligner(mode='global', substitution_matrix=blosum62, open_gap_score=-11, extend_gap_score=-1)

Local Alignment (Find Best Region)

aligner = PairwiseAligner(mode='local', match_score=2, mismatch_score=-1, open_gap_score=-10, extend_gap_score=-0.5)

Semiglobal (Overlap/Extension)

Allow free end gaps on query (useful for primer alignment)

aligner = PairwiseAligner(mode='global') aligner.query_left_open_gap_score = 0 aligner.query_left_extend_gap_score = 0 aligner.query_right_open_gap_score = 0 aligner.query_right_extend_gap_score = 0

Available Substitution Matrices

from Bio.Align import substitution_matrices print(substitution_matrices.load()) # List all available matrices

Common matrices

blosum62 = substitution_matrices.load('BLOSUM62') # General protein blosum80 = substitution_matrices.load('BLOSUM80') # Closely related proteins pam250 = substitution_matrices.load('PAM250') # Distantly related proteins

Working with SeqRecord Objects

from Bio import SeqIO

records = list(SeqIO.parse('sequences.fasta', 'fasta')) seq1, seq2 = records[0].seq, records[1].seq

aligner = PairwiseAligner(mode='global', match_score=1, mismatch_score=-1) alignments = aligner.align(seq1, seq2)

Iterating Over Multiple Alignments

Limit number of alignments returned (memory efficient)

aligner.max_alignments = 100

for i, alignment in enumerate(alignments): print(f'Alignment {i+1}: score={alignment.score}') if i >= 4: break

Substitution Matrix from Alignment

alignment = alignments[0] substitutions = alignment.substitutions

View as array (rows=target, cols=query)

print(substitutions)

Access specific substitution counts

substitutions['A', 'T'] gives count of A aligned to T

Export Alignment to Different Formats

alignment = alignments[0]

Various output formats

print(format(alignment, 'fasta')) # FASTA format print(format(alignment, 'clustal')) # Clustal format print(format(alignment, 'psl')) # PSL format (BLAT) print(format(alignment, 'sam')) # SAM format

Quick Reference: Scoring Parameters

Parameter Description Typical DNA Typical Protein

match_score

Score for identical bases 1-2 Use matrix

mismatch_score

Penalty for mismatches -1 to -3 Use matrix

open_gap_score

Cost to start a gap -5 to -15 -10 to -12

extend_gap_score

Cost per gap extension -0.5 to -2 -0.5 to -1

substitution_matrix

Scoring matrix N/A BLOSUM62

Common Errors

Error Cause Solution

OverflowError

Too many optimal alignments Set aligner.max_alignments

Low scores Wrong scoring scheme Use substitution matrix for proteins

No alignments in local mode Scores all negative Ensure match_score

0

Decision Tree: Choosing Alignment Mode

Need full-length comparison? ├── Yes → Use mode='global' │ └── Sequences similar length? │ ├── Yes → Standard global │ └── No → Consider semiglobal (free end gaps) └── No → Use mode='local' └── Find best matching regions only

Related Skills

• alignment-io - Save alignments to files in various formats

• msa-parsing - Work with multiple sequence alignments

• msa-statistics - Calculate identity, similarity metrics

• sequence-manipulation/motif-search - Pattern matching in sequences