Experiment 01Sequence SearchHeuristic Algorithms
Pairwise & Database Heuristic Sequence Search: NCBI BLAST (BLASTn, BLASTp), Bit Score Normalization & Karlin-Altschul E-Value Statistics
Duration: 4 Hours • Command-Line & Web BLAST
🎯 Experimental Aim & Learning Objectives
To perform local sequence similarity searches of unknown protein and nucleotide queries against NCBI reference databases using BLASTp and BLASTn; understand the 3-phase seed-extend heuristic algorithm; calibrate substitution matrices (BLOSUM62 vs PAM250); and mathematically evaluate statistical significance via raw alignment score ($S$), normalized bit score ($S’$), and the Karlin-Altschul Expectation value ($E$).
🔬 Theoretical Foundation & Statistical Equations
Exact dynamic programming algorithms (Smith-Waterman) are computationally prohibitive ($O(m times n)$) for querying massive genomic databases ($>10^{11}$ letters). The Basic Local Alignment Search Tool (BLAST) achieves orders-of-magnitude acceleration using a heuristic 3-phase strategy:
- Phase 1 (Seeding / Word Generation): The query sequence is broken into sliding words of length $W$ ($W=3$ amino acids for BLASTp, $W=11$ nucleotides for BLASTn). A neighborhood of similar words scoring $ge T$ when compared with the BLOSUM62 matrix is pre-computed.
- Phase 2 (Scanning & Two-Hit Triggering): The database is scanned for exact matches to the neighborhood words. Under the Two-Hit method, an extension is triggered only if two non-overlapping words occur on the same diagonal within window distance $A$.
- Phase 3 (Gapped Extension): The seed is extended bidirectionally without gaps, and then with a restricted dynamic programming band allowing insertions/deletions until the score drops by more than $X$ below the running maximum.
- Karlin-Altschul Statistical Model: The Expectation value ($E$) represents the number of distinct alignments with score $ge S$ expected to occur by chance alone in a random database search:
$$E = K cdot m cdot n cdot e^{-lambda S}$$
where $m$ is effective query length, $n$ is effective database size, and $K$ and $lambda$ are Karlin-Altschul extreme value distribution (Gumbel) scaling parameters. - Normalized Bit Score ($S’$): Converts raw score $S$ into an information-theoretic metric independent of scoring matrix or database scale:
$$S’ = frac{lambda S – ln(K)}{ln(2)} implies E = m cdot n cdot 2^{-S’}$$
A doubling of the database size requires an increase of 1 bit in $S’$ to maintain the same $E$-value threshold.
Figure 1.1: NCBI BLAST 3-phase seed-and-extend heuristic search architecture and Karlin-Altschul extreme value distribution statistics for raw alignment score ($S$), normalized bit score ($S’$), and expect value ($E$).
🧪 Computational Tools & Data Sets
- Platforms: NCBI Web BLAST Server (
blast.ncbi.nlm.nih.gov) and Local NCBI BLAST+ command-line suite (blastp,blastn,makeblastdbversion 2.14+). - Databases: RefSeq non-redundant protein database (
nr), Swiss-Prot/UniProtKB, and PDB sequence database. - Query Test Case: Unknown human hemoglobin subunit beta FASTA sequence (P68871) harboring the sickle-cell Glu6Val substitution.
📋 Step-by-Step Computational Protocol
- FASTA Acquisition: Retrieve the uncharacterized test FASTA sequence. Validate file formatting (single header line starting with
>, standard IUPAC single-letter amino acid code). - BLASTp Configuration: Open BLASTp. Select Database:
UniProtKB/Swiss-Prot(curated gold-standard). Set Algorithm parameters: Matrix =BLOSUM62, Gap Costs =Existence: 11, Extension: 1, Word size =6, Expect threshold =0.05. - Search Execution: Submit query. Inspect the Graphical Distribution of Hits, Hit Table, and pairwise High-Scoring Segment Pair (HSP) alignments.
- Matrix Sensitivity Comparison: Re-run the same query using
BLOSUM45(sensitive for distantly related homologs) andPAM30(strict for close orthologs). Observe how bit scores and identity percentages shift. - Command-Line Execution: Run equivalent search in bash:
blastp -query query.fasta -db swissprot -out results.tsv -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" -evalue 1e-5
📊 Representative Output & Statistical Evaluation
| Subject Accession | Description | Query Cover | Identity | Raw Score ($S$) | Bit Score ($S’$) | $E$-Value |
|---|---|---|---|---|---|---|
| P68871.2 | Hemoglobin subunit beta [Homo sapiens] | 100% | 99.3% | 774 | 301 bits | $4.0 times 10^{-105}$ |
| P02088.2 | Hemoglobin subunit beta-1 [Mus musculus] | 100% | 84.4% | 642 | 251 bits | $1.2 times 10^{-85}$ |
| P69905.2 | Hemoglobin subunit alpha [Homo sapiens] | 96% | 44.1% (Paralog) | 310 | 124 bits | $6.8 times 10^{-35}$ |
Biological Conclusion: The query is definitively identified as human hemoglobin beta chain with extreme statistical confidence ($E ll 10^{-100}$). The alpha subunit match demonstrates ancient gene duplication paralogy with significant sequence divergence.
{reading_nav}
