Sequence Alignment

Needleman-Wunsch performs global alignment, comparing two entire sequences end-to-end using dynamic programming with the recurrence F(i,j) = max{F(i-1,j-1)+s(xi,yj), F(i-1,j)+d, F(i,j-1)+d}. Smith-Waterman performs local alignment, finding the most similar subsequences by adding a zero option to the recurrence and starting traceback from the maximum score. Both have O(mn) time complexity but serve different purposes: global for homologous full-length sequences, local for finding conserved domains within divergent sequences.

BLAST uses a heuristic approach with four stages: (1) Seed word search identifies short exact matches (word_size, typically 11 for nucleotides, 3 for proteins). (2) Ungapped extension expands seeds in both directions while the score stays above threshold. (3) Gapped extension applies Smith-Waterman-like alignment to high-scoring segments. (4) Statistical evaluation calculates E-values using E = K*m*n*exp(-lambda*S). This pipeline is orders of magnitude faster than full Smith-Waterman while maintaining high sensitivity.

BLOSUM62 is directly derived from observed amino acid substitution frequencies in conserved blocks with 62% or more sequence identity, making it the standard for general protein database searches. PAM250 is constructed by extrapolating from the PAM1 matrix (1% substitution probability) raised to the 250th power, using an evolutionary model. PAM250 has similar sensitivity to BLOSUM45 and is better suited for detecting distant evolutionary relationships. BLOSUM matrices are empirical; PAM matrices are theoretical.

ClustalW/Omega uses a progressive method building alignments along a guide tree from pairwise distances. MUSCLE uses iterative refinement to improve an initial progressive alignment, offering better accuracy than ClustalW. MAFFT provides multiple strategies: --auto for automatic selection, --localpair (L-INS-i) for highest accuracy with fewer sequences, and default mode for speed. In general, speed ranking is ClustalW < MUSCLE < MAFFT, and MAFFT L-INS-i typically provides the best accuracy for small to medium datasets.

The E-value represents the expected number of alignments with that score occurring by chance in a database of the given size. E-values below 1e-50 indicate definite homology, below 1e-5 indicate statistically significant similarity, and above 1 suggest possible random matches. The formula is E = K*m*n*exp(-lambda*S), where m is query length, n is database size, and S is the raw alignment score. Smaller E-values indicate more biologically meaningful alignments.

Linear gap penalty applies a constant cost per gap position: W(k) = -d*k. Affine gap penalty distinguishes between opening a new gap and extending an existing one: W(k) = -d - e*(k-1), where d is the gap open penalty and e is the gap extension penalty. Typical values with BLOSUM62 are d=10, e=0.5. Affine penalties are biologically more realistic because insertions and deletions in evolution tend to occur as contiguous blocks rather than isolated single-residue events.

FASTA format uses > headers followed by sequence lines (most universal). CLUSTAL format shows aligned sequences with conservation symbols (*, :, .) below columns. Stockholm format (used by Pfam/Rfam) supports rich annotations including secondary structure consensus. Phylip format uses fixed 10-character sequence names with a header line specifying sequence count and alignment length. Most alignment tools can read and write multiple formats, but Stockholm is preferred for databases and FASTA for general use.

Sum-of-Pairs (SP) score sums all pairwise substitution scores at each column, with higher SP indicating better alignment. Total Column (TC) score measures the fraction of columns where all sequences agree. Shannon entropy (H = -sum(pi*log2(pi))) quantifies conservation at each position, where H=0 means fully conserved. Dot plots provide visual assessment of pairwise similarity. Tools like Jalview display conservation histograms and consensus sequences for interactive quality evaluation.