Day 7: Pairwise Sequence Alignment, Scoring Matrices & Heuristic Database Searching (BLAST)

1. Dynamic Programming Foundations in Sequence Alignment

Sequence alignment is the fundamental computational method for establishing homology, identifying conserved active sites, and tracing molecular evolution. Alignment algorithms fall into two broad algorithmic classes:

2. Amino Acid Substitution Scoring Matrices: PAM vs. BLOSUM

While nucleotide alignment uses simple match (+1) and mismatch (-2) scores, protein alignment relies on empirical log-odds substitution matrices that reflect chemical similarity, mutational frequencies, and genetic code architecture:

$$\text{Score}(a, b) = 10 \cdot \log_{10} \left( \frac{q_{ab}}{p_a \cdot p_b} \right)$$

PAM (Point Accepted Mutation) Matrices

• Constructed by Margaret Dayhoff based on 71 families of very closely related proteins (> 85% sequence identity).
• Higher PAM numbers correspond to greater evolutionary distance: PAM1 corresponds to 1 point mutation per 100 amino acids (1% divergence); PAM250 corresponds to 250 mutations per 100 amino acids (~20% sequence identity remaining).

BLOSUM (Blocks Substitution Matrix) Matrices

• Constructed by Steven & Jorja Henikoff derived from ungapped local alignments of conserved protein blocks in the BLOCKS database.
• Lower BLOSUM numbers correspond to greater evolutionary distance:
  • BLOSUM62: Default matrix for general-purpose searches; sequences with $\ge 62\%$ identity are clustered. Ideal for detecting moderately divergent homologues (30–60% identity).
  • BLOSUM80: Clustered at 80% identity. Stricter matrix optimized for closely related homologues.
  • BLOSUM45: Clustered at 45% identity. Permissive matrix optimized for highly divergent, distant homologues (< 30% identity).

3. Heuristic Searching & Karlin-Altschul BLAST Statistics

While Smith-Waterman is mathematically optimal ($O(mn)$ time complexity), database searches across hundreds of millions of sequences require heuristic speedups. The Basic Local Alignment Search Tool (BLAST) achieves orders-of-magnitude speedups by indexing short "words" (seeds: $W=11$ for blastn, $W=3$ for blastp), extending matches via two-hit ungapped alignments, and performing gapped band-limited extensions.

Karlin-Altschul E-Value Formula

The statistical significance of an alignment score $S$ is quantified by the Expectation Value ($E$-value)—the number of distinct alignments with score $\ge S$ expected to occur purely by chance when searching a database of size $N$ with a query of length $m$:

$$E = K \cdot m \cdot N \cdot e^{-\lambda S}$$

where $K$ and $\lambda$ are natural scale parameters derived from the scoring system and background amino acid frequencies. The normalized Bit Score ($S'$) eliminates database and scoring system dependencies:

$$S' = \frac{\lambda S - \ln K}{\ln 2} \implies E = m \cdot N \cdot 2^{-S'}$$

4. Production Bash & Python Script: Local BLAST+ Execution

#!/usr/bin/env bash
# Day 7 Practical: Local BLAST+ Database Construction & High-Throughput Search

# 1. Build a custom local BLAST database from a FASTA proteome
makeblastdb -in uniref50_sample.fasta -dbtype prot -out custom_proteome_db -title "ReferenceProteome"

# 2. Run blastp with standard tabular output format (-outfmt 6)
# Format fields: query id, subject id, % identity, alignment length, mismatches, gap opens,
# q. start, q. end, s. start, s. end, evalue, bit score
blastp -query query_protein.fasta        -db custom_proteome_db        -evalue 1e-5        -max_target_seqs 10        -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore"        -out blastp_results.tsv

echo "BLASTp search complete. Results written to blastp_results.tsv"