🧠 BLAST Fundamentals

Q: What is BLAST?

BLAST (Basic Local Alignment Search Tool) is a heuristic algorithm for comparing biological sequences (DNA, RNA, or protein) against sequence databases. It finds regions of local similarity by using a seed-and-extend strategy, trading some sensitivity for dramatic speed improvements over exact algorithms like Smith-Waterman.

Q: What does E-value mean in BLAST?

The E-value (Expectation value) is the number of hits with a score at least as good as the observed score that would be expected by chance when searching a database of a particular size. Smaller E-values indicate more significant matches. E-values below 1e-5 are typically considered significant; below 1e-50 are highly significant.

Q: What is the difference between BLAST and Smith-Waterman?

BLAST is a heuristic algorithm that uses seed-and-extend to find candidate alignments quickly — it's 100-1000x faster than Smith-Waterman. Smith-Waterman is exact and guaranteed optimal but much slower. BLAST may miss some remote homologs that Smith-Waterman would find. Use BLAST for database searches; use Smith-Waterman for confirmation and high-sensitivity alignment.

Q: What are the 5 BLAST variants?

The five main BLAST programs are: (1) blastn — nucleotide query vs nucleotide database; (2) blastp — protein query vs protein database; (3) blastx — translated nucleotide query vs protein database; (4) tblastn — protein query vs translated nucleotide database; (5) tblastx — translated nucleotide query vs translated nucleotide database.

Q: What is a bit score in BLAST?

The bit score (S') is a normalized alignment score that is independent of database size. It is derived from the raw alignment score (S) using the statistical parameters of the scoring system: S' = (λS - ln K) / ln 2. Higher bit scores indicate better alignments. Bit scores can be compared across different database searches, unlike raw scores.

🎯 What is BLAST?

BLAST (Basic Local Alignment Search Tool) is a heuristic algorithm for comparing biological sequences (DNA, RNA, or protein) against sequence databases. Developed by Stephen Altschul and colleagues at NCBI in 1990, it uses a seed-and-extend strategy to find regions of local similarity, trading some sensitivity for 100-1000x speed improvements over exact algorithms like Smith-Waterman. It's the most widely used tool in bioinformatics.

Why BLAST Matters

Speed: 100-1000x faster than SW. Sensitive: Finds homologs. Statistical: E-values. Variants: 5 programs. Ubiquitous: Every biologist uses it. Foundation: Basis for DIAMOND, MMseqs2.

📌 Heuristic📌 Seed-and-Extend📌 E-value📌 Database Search

# BLAST — the big picture
# Basic Local Alignment Search Tool.

# CONCEPTUAL OVERVIEW
print("""
=== BLAST ===

What it is:
  - Heuristic sequence search
  - Local alignment
  - Database search
  - Fast (100-1000x vs SW)
  - Statistical significance

Created by:
  - Stephen Altschul et al. (1990)
  - NCBI

Key concepts:
  - Seed: Short exact match
  - Extend: Grow alignment
  - Score: Alignment score
  - Bit score: Normalized
  - E-value: Significance
  - HSP: High-scoring Segment Pair

Variants (5):
  - blastn: Nucleotide vs nucleotide
  - blastp: Protein vs protein
  - blastx: Translated nucleotide vs protein
  - tblastn: Protein vs translated nucleotide
  - tblastx: Translated vs translated

What it's NOT:
  - Not exact (heuristic)
  - Not always optimal
  - Not for structure
  - Not for phylogeny (alone)

What it IS:
  - Database search
  - Homology detection
  - Function prediction
  - Foundation of genomics
""")

# WHY IT MATTERS
print("=== Why It Matters ===")
print("""
Biological questions:
  - What is this gene?
  - What does this protein do?
  - Is this sequence novel?
  - Are these sequences related?

BLAST answers:
  - Finds similar sequences
  - Predicts function
  - Identifies homologs
  - Estimates significance
  - Fast and scalable
""")

# REAL-WORLD SCALE
print("=== Real-World Scale ===")
print("""
NCBI BLAST:
  - Millions of searches per day
  - 100+ databases
  - Petabytes of sequence data

Databases:
  - nr (non-redundant)
  - nt (nucleotide)
  - RefSeq
  - SwissProt
  - PDB
  - UniProt

Usage:
  - Every molecular biology lab
  - Every genomics project
  - Every bioinformatics tool
  - Foundation of biology
""")

🔀 BLAST Variants

🎯 Five Programs

BLAST has five main variants for different combinations of query and database types: nucleotide or protein, direct or translated.

The 5 Variants

blastn: DNA vs DNA. blastp: Protein vs protein. blastx: DNA vs protein. tblastn: Protein vs DNA. tblastx: DNA vs DNA (translated).

🔀 blastn🔀 blastp🔀 blastx🔀 tblastn🔀 tblastx

# BLAST Variants
print("=== 1. blastn ===")
print("""
Nucleotide query vs Nucleotide database

Query: DNA/RNA sequence
Database: DNA/RNA sequences
Comparison: Direct

Use cases:
  - Identify a gene
  - Find similar genomes
  - Primer specificity
  - SNP analysis

Example:
  blastn -query gene.fasta -db nt -out results.txt

Speed: Fast
Sensitivity: Medium
""")

print("\n=== 2. blastp ===")
print("""
Protein query vs Protein database

Query: Protein sequence
Database: Protein sequences
Comparison: Direct

Use cases:
  - Protein function
  - Homology
  - Domain detection
  - Protein family

Example:
  blastp -query protein.fasta -db nr -out results.txt

Speed: Medium
Sensitivity: High
""")

print("\n=== 3. blastx ===")
print("""
Translated Nucleotide query vs Protein database

Query: DNA/RNA (translated in 6 frames)
Database: Protein sequences
Comparison: Translated query

Use cases:
  - Annotate new genes
  - EST analysis
  - Metagenomics
  - Find coding regions

Example:
  blastx -query dna.fasta -db nr -out results.txt

Speed: Slower (6 frames)
Sensitivity: High
""")

print("\n=== 4. tblastn ===")
print("""
Protein query vs Translated Nucleotide database

Query: Protein sequence
Database: DNA (translated in 6 frames)
Comparison: Translated database

Use cases:
  - Find genes in genomes
  - Identify ORFs
  - Comparative genomics
  - Homology in DNA

Example:
  tblastn -query protein.fasta -db nt -out results.txt

Speed: Slower (6 frames)
Sensitivity: High
""")

print("\n=== 5. tblastx ===")
print("""
Translated Nucleotide query vs Translated Nucleotide database

Query: DNA (translated in 6 frames)
Database: DNA (translated in 6 frames)
Comparison: Translated both

Use cases:
  - Distant homology
  - Novel genes
  - Metagenomics
  - Deep evolution

Example:
  tblastx -query dna.fasta -db nt -out results.txt

Speed: Slowest (36 comparisons)
Sensitivity: Highest
""")

print("\n=== Summary Table ===")
print(f"{'Variant':<12} {'Query':<15} {'Database':<15} {'Frames':<10}")
print("-" * 60)
print(f"{'blastn':<12} {'DNA':<15} {'DNA':<15} {'1':<10}")
print(f"{'blastp':<12} {'Protein':<15} {'Protein':<15} {'1':<10}")
print(f"{'blastx':<12} {'DNA':<15} {'Protein':<15} {'6':<10}")
print(f"{'tblastn':<12} {'Protein':<15} {'DNA':<15} {'6':<10}")
print(f"{'tblastx':<12} {'DNA':<15} {'DNA':<15} {'36':<10}")

print("\n=== Which Variant to Use? ===")
print("""
Use blastn when:
  - Both sequences are DNA
  - Looking for exact or near-exact matches
  - Same species or close species

Use blastp when:
  - Both sequences are proteins
  - Most sensitive for protein homology
  - Annotating proteins

Use blastx when:
  - Query is DNA, want to find protein homologs
  - Annotating new genomic DNA
  - EST/mRNA analysis

Use tblastn when:
  - Query is protein, search DNA database
  - Finding genes in genome
  - Searching draft genomes

Use tblastx when:
  - Both are DNA, want max sensitivity
  - Distant homology
  - Novel genes
  - Slow but thorough
""")

⚙️ The Algorithm

🎯 How BLAST Works

BLAST uses a three-step heuristic: (1) find short seeds (words), (2) extend them into high-scoring pairs (HSPs), (3) evaluate statistical significance. This is much faster than Smith-Waterman but may miss some alignments.

The Algorithm

1. Seeding: Find short matches. 2. Extension: Grow seeds. 3. Evaluation: Score and filter. Word size: 11 (DNA), 3 (protein). T threshold: Score threshold. HSP: High-scoring pair.

⚙️ Seed⚙️ Extend⚙️ Evaluate⚙️ HSP

# The Algorithm
print("=== Three Steps ===")
print("""
1. SEEDING
   - Break query into k-mers (words)
   - k = 11 (nucleotide), 3 (protein)
   - Find exact matches in database
   - Also near-matches (using neighborhood)

2. EXTENSION
   - Extend seeds both directions
   - Use Smith-Waterman-like DP
   - Stop when score drops below threshold
   - Generate HSPs (High-scoring Segment Pairs)

3. EVALUATION
   - Compute bit score
   - Compute E-value
   - Rank by significance
   - Report top hits
""")

print("\n=== Word Size ===")
print("""
Nucleotide (blastn):
  - Word size: 11 (default)
  - 4^11 = 4M possible words
  - Random match: rare
  - Sensitive

Protein (blastp):
  - Word size: 3 (default)
  - 20^3 = 8000 words
  - Random match: more common
  - Neighborhood: similar amino acids

Trade-off:
  - Larger word: faster, less sensitive
  - Smaller word: slower, more sensitive

Tuning:
  - blastn: 7-11 words
  - blastp: 2-7 words
  - Use smaller for distant homology
""")

print("\n=== Seeding Example ===")
print("""
Query: ATGCGTACGTTAG

Words (k=4):
  ATGC
  TGCG
  GCGT
  CGTA
  GTAC
  TACG
  ACGT
  CGTT
  GTTA
  TTAG

Database:
  ATGCGTAC...
  ^^^^
  Match found!
  
Then extend in both directions.
""")

print("\n=== Extension ===")
print("""
Seeds extended with DP:
  - Use Smith-Waterman-like
  - Stop when score drops
  - X-dropoff: stop if score drops by X
  - Default X = 15-30

Example:
  Seed: ATGC
  Extend right: ATGCGTACG...
  Extend left: ...ATGC
  
  Continue while score improves.
""")

print("\n=== HSPs ===")
print("""
HSP: High-scoring Segment Pair

Definition:
  - Ungapped alignment
  - Two sequences
  - Score above threshold

Multiple HSPs:
  - Same query, different DB
  - Same query, same DB, different regions
  
Combine into final alignment.

Statistics:
  - Bit score
  - E-value
""")

print("\n=== Detailed Algorithm ===")
print("""
1. INDEX DATABASE
   - Precompute k-mer index
   - Once per database
   - Fast queries

2. FIND SEEDS
   - For each query k-mer
   - Look up in index
   - Get positions

3. EXTEND SEEDS
   - Both directions
   - Score with matrix
   - X-dropoff
   - Generate HSPs

4. FILTER HSPs
   - Score threshold (S)
   - Statistical significance
   - Remove low-quality

5. COMBINE HSPs
   - Same query-database pair
   - Gapped alignment
   - Final HSP

6. REPORT
   - Sort by E-value
   - Top hits
   - Alignments
""")

print("\n=== Complexity ===")
print("""
Time: O(m + n × k) where k is word count
  - Much faster than SW O(mn)
  - For database search: O(m × D) where D = database size
  - Practically: seconds to minutes

Space: O(n) for index
  - Precomputed
  - Shared across queries

Speedup vs SW:
  - 100-1000x
  - Depends on similarity
  - More similar = faster
""")

🌱 Seed-and-Extend

🎯 The Core Heuristic

Seed-and-extend is BLAST's key optimization: find short exact matches (seeds), then extend them into longer alignments. This avoids the O(mn) cost of full DP.

Seed-and-Extend

Seeds: Short exact words. Neighborhood: Similar words. Extension: Both directions. X-dropoff: Stop criterion. HSPs: Output. Tandem: Combine.

🌱 Seeds🌱 Neighborhood🌱 Extension🌱 X-dropoff

# Seed-and-Extend
print("=== The Idea ===")
print("""
Observation:
  - Homologous sequences share short exact matches
  - These seeds can be found quickly
  - Extend to find longer alignments

Steps:
  1. Build index of database words
  2. For each query word, find matches
  3. Extend matches into HSPs
  4. Combine and report
""")

print("\n=== Exact Seeds ===")
print("""
For nucleotides:
  - Word size 11 (default)
  - Exact match required
  - 4^11 = 4M possible words
  - Random match probability: 1/4M

For proteins:
  - Word size 3 (default)
  - Neighborhood of similar words
  - BLOSUM62 defines similarity
  - e.g., "LIV" neighbors: "LIV", "LVI", etc.
""")

print("\n=== Neighborhood Words ===")
print("""
Protein neighborhood:
  - Exact word: LIV
  - Neighborhood: similar words
  - Threshold T: score cutoff

Example (T=11, BLOSUM62):
  LIV (12)
  LVI (11)
  ILV (12)
  IVL (11)
  VLI (11)
  VIL (11)
  ... plus similar
  
Total: ~50 words

Why?
  - Conserved substitutions
  - Biological similarity
  - More sensitive
""")

print("\n=== Extension ===")
print("""
Extend seeds:
  - Both directions (5' and 3')
  - Use scoring matrix
  - Add matches/mismatches/gaps
  - Stop when score drops

X-dropoff:
  - X = 15-30 typically
  - If score drops by X, stop
  - Prevents extending into noise
  - Fast termination
""")

print("\n=== Example ===")
print("""
Query: ATGCGTACGTTAG
DB:    ATGCGTAGTTAG

Step 1: Find seeds (k=4)
  ATGC matches position 0
  TGCG matches position 1
  ... etc.

Step 2: Extend
  From ATGC, extend right:
    A: score +1 (match)
    T: +1
    G: +1
    C: +1
    G: +1
    T: +1
    A: +1
    C: +1
    G: +1
    T: +1
    T: -1 (mismatch: C vs G? wait...)
    
  Let me redo:
  Query: ATGCGTACGTTAG
  DB:    ATGCGTAGTTAG
         ^^^^^^^^ 
         Matches: 12
         Mismatch: 1
         
Step 3: HSP score = 12 - 1 = 11 (with match=1, mismatch=-1)

Step 4: E-value computed
""")

print("\n=== Two-Hit Heuristic ===")
print("""
BLAST uses two-hit heuristic:
  - Find two non-overlapping HSPs on same diagonal
  - Combine them
  - More efficient than single-hit
  - Reduces false positives

Benefit:
  - 3x faster than single-hit
  - Same sensitivity
  - Standard in BLAST
""")

print("\n=== Index Structures ===")
print("""
Suffix array:
  - Sorted suffixes
  - Fast lookup
  - Memory efficient

FM-index:
  - Burrows-Wheeler transform
  - Compressed
  - Used in BWA, Bowtie

Hash table:
  - k-mer → positions
  - Used in BLAST
  - Fast but memory heavy

Choice:
  - BLAST: hash table
  - Modern aligners: FM-index
""")

print("\n=== Optimization ===")
print("""
1. Word size
   - Larger = faster, less sensitive
   - Smaller = slower, more sensitive

2. Neighborhood size
   - Larger = more sensitive
   - Smaller = faster

3. X-dropoff
   - Larger = more sensitive
   - Smaller = faster

4. Two-hit
   - Faster
   - Standard

5. Database indexing
   - Precompute
   - Reuse
""")

📊 Statistics & E-value

🎯 Measuring Significance

BLAST uses Karlin-Altschul statistics to evaluate the significance of alignments. The two key measures are the bit score (normalized score) and the E-value (expected number of random hits).

Statistics

Bit score: Normalized score. E-value: Expected hits. Karlin-Altschul: Theory. λ, K: Parameters. P-value: Probability. Smaller: Better.

📊 Bit Score📊 E-value📊 Karlin-Altschul📊 Significance

# Statistics & E-value
print("=== Raw Score ===")
print("""
Raw score (S):
  - Sum of substitution scores
  - Sum of gap penalties
  - Not normalized
  - Depends on matrix

Example:
  Match: +1
  Mismatch: -1
  Gap: -1
  
  Alignment:
    ATGCGT
    ATGCGT
    Score: 6

Not comparable across matrices.
""")

print("\n=== Bit Score ===")
print("""
Bit score (S'):
  - Normalized score
  - Independent of matrix
  - Comparable across searches

Formula:
  S' = (λ × S - ln K) / ln 2

Where:
  λ = scale parameter
  K = scale parameter
  S = raw score

Properties:
  - Always positive
  - Higher = better
  - Doubling score = 1 more bit
""")

print("\n=== E-value ===")
print("""
E-value:
  - Expected number of random hits
  - With score ≥ S in database
  
Formula:
  E = K × m × n × e^(-λS)

Where:
  m = query length
  n = database length
  S = raw score

Interpretation:
  E = 1e-5: Very significant
  E = 1e-10: Highly significant
  E = 0.01: Marginal
  E = 1: Random
  E = 10: Not significant
""")

print("\n=== E-value Examples ===")
print("""
E = 1e-100: Extremely significant
  - Essentially identity
  - Same gene

E = 1e-50: Very significant
  - Clear homolog
  - Same protein family

E = 1e-20: Significant
  - Likely homolog
  - Similar function

E = 1e-5: Marginally significant
  - Possible homolog
  - Worth investigating

E = 0.01: Marginal
  - Maybe homologous
  - Need more evidence

E = 1: Random
  - Not significant

E = 10: Definitely random
""")

print("\n=== Bit Score vs E-value ===")
print("""
Relationship:
  E = m × n × 2^(-S')

Where:
  m = query length
  n = database length
  S' = bit score

At E = 1:
  S' = log2(m × n)

Example:
  m = 100, n = 1e9
  S' = log2(1e11) ≈ 37 bits

Rule of thumb:
  E < 1e-5: significant
  E < 1e-10: highly significant
  E < 1e-50: essentially identical
""")

print("\n=== P-value ===")
print("""
P-value:
  - Probability of observing such a hit by chance
  - E = -ln(1-P) → P = 1 - e^(-E)
  
For small E:
  P ≈ E

Interpretation:
  Same as E-value
  - P < 1e-5 significant
""")

print("\n=== Karlin-Altschul Theory ===")
print("""
Requirements:
  1. Random sequences
  2. Scoring matrix with negative expected score
  3. Ungapped alignment

Results:
  - Scores follow extreme value distribution
  - E-value formula derived
  - Parameters λ, K depend on matrix

For gapped:
  - More complex
  - Empirical parameters
  - Still works

Important:
  - Assumes random sequences
  - Biological sequences not random
  - E-values are approximations
""")

print("\n=== Significance Guidelines ===")
print(f"{'E-value':<15} {'Interpretation':<40}")
print("-" * 55)
print(f"{'< 1e-100':<15} {'Essentially identical':<40}")
print(f"{'1e-50 to 1e-100':<15} {'Same protein family':<40}")
print(f"{'1e-20 to 1e-50':<15} {'Clear homolog':<40}")
print(f"{'1e-5 to 1e-20':<15} {'Likely homolog':<40}")
print(f"{'1e-3 to 1e-5':<15} {'Marginal':<40}")
print(f"{'> 1e-3':<15} {'Not significant':<40}")

print("\n=== Practical Tips ===")
print("""
1. Use E-value threshold
   - Default: 10
   - Strict: 1e-5
   - Very strict: 1e-10

2. Consider database size
   - Larger database → higher E
   - Adjust expectations

3. Bit score is portable
   - Compare across searches
   - E-value depends on DB size

4. Don't ignore marginal hits
   - E = 0.01 might be real
   - Look at alignment

5. Query coverage matters
   - Partial alignments
   - Domain hits

6. Multiple HSPs
   - Same hit, different regions
   - Combine for full picture
""")

📐 Scoring Matrices

🎯 Scoring Similarity

BLAST uses substitution matrices to score matches and mismatches. For proteins, BLOSUM62 is standard. For nucleotides, simpler scoring is used.

Scoring

BLOSUM62: Protein default. PAM: Alternative. Match/mismatch: DNA. Gap: Open + extend. Context: Depends on distance. Trade-off: Sensitivity vs speed.

📐 BLOSUM62📐 PAM📐 Match/Mismatch📐 Gap

# Scoring Matrices
print("=== Protein Matrices ===")
print("""
BLOSUM (BLOcks SUbstitution Matrix):
  - BLOSUM45: Distant (< 45% identity)
  - BLOSUM50: Distant
  - BLOSUM62: General (default)
  - BLOSUM80: Close (> 80%)
  - BLOSUM90: Very close

PAM (Point Accepted Mutation):
  - PAM30: Distant
  - PAM70: Distant
  - PAM250: Very distant (default for some)

Recommendation:
  - BLOSUM62 for most cases
  - Lower BLOSUM for distant homologs
  - Higher BLOSUM for close homologs
""")

print("\n=== BLOSUM62 (partial) ===")
print("""
      A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
A     4  -1  -2  -2   0  -1  -1   0  -2  -1  -1  -1  -1  -2  -1   1   0  -3  -2   0
R    -1   5   0  -2  -3   1   0  -2   0  -3  -2   2  -1  -3  -2  -1  -1  -3  -2  -3
N    -2   0   6   1  -3   0   0   0   1  -3  -3   0  -2  -3  -2   1   0  -4  -2  -3
D    -2  -2   1   6  -3   0   2  -1  -1  -3  -4  -1  -3  -3  -1   0  -1  -4  -3  -3
C     0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
...
""")

print("\n=== DNA Scoring ===")
print("""
Simple:
  Match: +1
  Mismatch: -1
  Gap: -2

Standard (BLAST):
  Match: +2
  Mismatch: -3
  Gap open: -5
  Gap extend: -2

Transition/transversion:
  Transition (A↔G, C↔T): +1
  Transversion: -3
""")

print("\n=== Choosing a Matrix ===")
print(f"{'Scenario':<30} {'Matrix':<20} {'Gap':<20}")
print("-" * 70)
print(f"{'General protein':<30} {'BLOSUM62':<20} {'11+1':<20}")
print(f"{'Close protein':<30} {'BLOSUM80':<20} {'11+1':<20}")
print(f"{'Distant protein':<30} {'BLOSUM45':<20} {'11+1':<20}")
print(f"{'Very distant':<30} {'BLOSUM45':<20} {'14+2':<20}")
print(f"{'DNA standard':<30} {'match+2/-3':<20} {'5+2':<20}")
print(f"{'DNA distant':<30} {'match+1/-1':<20} {'2+2':<20}")

print("\n=== Example Matrix ===")
print("""
# BLOSUM62 (Python)
BLOSUM62 = {
    'A': {'A': 4, 'R': -1, 'N': -2, 'D': -2, 'C': 0, 'Q': -1, 'E': -1, 'G': 0, 
          'H': -2, 'I': -1, 'L': -1, 'K': -1, 'M': -1, 'F': -2, 'P': -1, 'S': 1, 
          'T': 0, 'W': -3, 'Y': -2, 'V': 0},
    'R': {'A': -1, 'R': 5, 'N': 0, 'D': -2, 'C': -3, 'Q': 1, 'E': 0, 'G': -2,
          'H': 0, 'I': -3, 'L': -2, 'K': 2, 'M': -1, 'F': -3, 'P': -2, 'S': -1,
          'T': -1, 'W': -3, 'Y': -2, 'V': -3},
    # ... (all 20 amino acids)
}

def score(a, b, matrix=BLOSUM62):
    return matrix[a][b]
""")

print("\n=== Gap Penalties ===")
print("""
Linear:
  gap = -g × length
  
Affine:
  gap = -(open + extend × (length - 1))

BLAST defaults:
  Protein:
    - BLOSUM62: 11 + 1
    - BLOSUM45: 14 + 2
    - BLOSUM80: 10 + 1
  DNA:
    - match +2/-3: 5 + 2
    - match +1/-1: 2 + 2
""")

➖ Gap Penalties

🎯 Modeling Indels

BLAST uses affine gap penalties to model insertions and deletions realistically: opening a gap is expensive, extending it is cheap.

Gap Penalties

Affine: Open + extend. Open: Cost to start. Extend: Cost per residue. Default: 11+1 (protein). DNA: 5+2. Biological: Realistic.

➖ Affine➖ Open➖ Extend➖ Biological

# Gap Penalties
print("=== Why Affine? ===")
print("""
Biological reality:
  - One long indel > many short ones
  - DNA replication slippage
  - Recombination
  
Affine penalty:
  gap = -(open + extend × (length - 1))

Example (open=11, extend=1):
  Length 1: -11
  Length 2: -12
  Length 3: -13
  Length 10: -20
  
Compare to linear (gap=5):
  Length 1: -5
  Length 10: -50
""")

print("\n=== BLAST Default Gap Penalties ===")
print(f"{'Matrix':<20} {'Open':<10} {'Extend':<10}")
print("-" * 45)
print(f"{'BLOSUM62':<20} {'11':<10} {'1':<10}")
print(f"{'BLOSUM45':<20} {'14':<10} {'2':<10}")
print(f"{'BLOSUM80':<20} {'10':<10} {'1':<10}")
print(f"{'PAM30':<20} {'9':<10} {'1':<10}")
print(f"{'PAM70':<20} {'10':<10} {'1':<10}")
print(f"{'DNA +2/-3':<20} {'5':<10} {'2':<10}")
print(f"{'DNA +1/-1':<20} {'2':<10} {'2':<10}")

print("\n=== Effect of Gap Penalties ===")
print("""
Low gap penalty:
  - More gaps
  - Shorter alignments
  - More sensitivity
  - More false positives

High gap penalty:
  - Fewer gaps
  - Longer alignments
  - Less sensitivity
  - Fewer false positives

Balance:
  - Match score > gap open
  - Extend < open
  - Tuned for biology
""")

print("\n=== Gotoh's Algorithm ===")
print("""
Gotoh's algorithm:
  - Implements affine gaps
  - O(mn) time (same as linear)
  - Three matrices:
    M: match/mismatch
    Ix: gap in x
    Iy: gap in y

Used by:
  - BLAST (for gapped alignment)
  - Most aligners
""")

print("\n=== Practical Example ===")
print("""
Sequences:
  A: ATGCGTACGTTAG
  B: ATGCGT---GTTAG
              (3bp deletion)

Score with affine (open=5, extend=1):
  Matches: 12 × 2 = 24
  Gap: -(5 + 1×2) = -7
  Total: 17

Score with linear (gap=5):
  Matches: 24
  Gap: -5 × 3 = -15
  Total: 9

Affine rewards the single long gap.
""")

print("\n=== When to Adjust ===")
print("""
Increase gap penalty when:
  - Sequences very similar
  - Expect few indels
  - Reduce false positives

Decrease gap penalty when:
  - Distant homologs
  - Expect many indels
  - Increase sensitivity

Default is usually best.
""")

📄 Output Format

🎯 Reading BLAST Results

BLAST output has a standard format with three sections: header (program info), hits list (ranked matches), and alignments (detailed HSPs).

Output

Header: Program, database, query. Hits: Ranked list. Alignments: Detailed. Score: Bit score. E-value: Significance. Identity: Percent.

📄 Header📄 Hits📄 Alignments📄 Score

# Output Format
print("=== Standard BLAST Output ===")
print("""
BLASTN 2.13.0+
Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb Miller (2000)

Query= query_sequence
Length=100

Database: nt
           100,000 sequences; 100,000,000 total letters

Results:
                                                                     Score    E
Sequences producing significant alignments:                        (bits) Value

ref|NC_000001.11| Homo sapiens chromosome 1                        200    1e-50
ref|NC_000002.12| Homo sapiens chromosome 2                        180    1e-45
ref|NC_000003.12| Homo sapiens chromosome 3                        150    1e-35

>ref|NC_000001.11| Homo sapiens chromosome 1
Length=248956422

 Score = 200 bits (108),  Expect = 1e-50
 Identities = 100/100 (100%), Gaps = 0/100 (0%)
 Strand=Plus/Minus

Query  1    ATGCGTACGTTAG...  60
             |||||||||||||
Sbjct  100  ATGCGTACGTTAG...  41
""")

print("\n=== Sections ===")
print("""
1. HEADER
   - Program version
   - Reference
   - Query name
   - Query length
   - Database name
   - Database size

2. HITS LIST
   - Ranked by E-value
   - Score (bits)
   - E-value
   - Accession
   - Description

3. ALIGNMENTS
   - Detailed HSPs
   - Query sequence
   - Match line
   - Subject sequence
   - Strand
   - Position
""")

print("\n=== Key Fields ===")
print("""
Score (bits):
  - Normalized score
  - Higher = better
  - Comparable across searches

Expect (E-value):
  - Expected random hits
  - Lower = better
  - Threshold for significance

Identities:
  - Exact matches / total
  - Percent identity

Gaps:
  - Gap positions / total
  - Percent gaps

Strand:
  - Plus/Plus or Plus/Minus
  - For DNA

Frame:
  - Reading frame
  - For translated
""")

print("\n=== Output Formats ===")
print("""
1. Pairwise (default)
   - Human-readable
   - Alignments

2. Tabular (outfmt 6)
   - Tab-separated
   - Machine-readable
   - Standard

3. XML (outfmt 5)
   - Structured
   - Parseable
   - Verbose

4. JSON (outfmt 15)
   - Modern
   - Parseable

5. CSV (outfmt 10)
   - Comma-separated
""")

print("\n=== Tabular Format (outfmt 6) ===")
print("""
Default columns:
  1. qseqid: Query seq ID
  2. sseqid: Subject seq ID
  3. pident: % identity
  4. length: Alignment length
  5. mismatch: Mismatches
  6. gapopen: Gap openings
  7. qstart: Query start
  8. qend: Query end
  9. sstart: Subject start
  10. send: Subject end
  11. evalue: E-value
  12. bitscore: Bit score

Custom:
  -outfmt "6 qseqid sseqid pident length evalue bitscore"
""")

print("\n=== Example ===")
print("""
query1  ref|NC_000001.11|  100.000  100  0  0  1  100  100  1  1e-50  200
query1  ref|NC_000002.12|  95.000   100  5  0  1  100  100  1  1e-45  180
query1  ref|NC_000003.12|  90.000   100  10 0  1  100  100  1  1e-35  150
""")

print("\n=== Reading Output ===")
print("""
1. SORT BY E-VALUE
   - Most significant first

2. CHECK IDENTITY
   - High identity = close homolog
   - Low identity = distant

3. CHECK COVERAGE
   - Full-length = likely homolog
   - Partial = domain

4. CHECK STRAND
   - Plus/Minus for reverse complement

5. MULTIPLE HITS
   - Same subject = multiple HSPs
   - Combine for full alignment
""")

🐍 Python Integration

🎯 Running BLAST from Python

You can run BLAST from Python using Biopython (local or remote) or parse existing BLAST output. Here's how.

Python Integration

Biopython: NCBIWWW, NcbiblastnCommandline. Remote: NCBI BLAST API. Local: blastn binary. Parse: NCBIXML, SearchIO. Automate: Batch, pipeline.

🐍 Biopython🐍 NCBIWWW🐍 Parse🐍 Automate

# Python Integration
print("=== Remote BLAST (NCBI) ===")
print("""
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML
from Bio import SeqIO

# Read query
record = SeqIO.read("query.fasta", "fasta")

# Run BLAST
result_handle = NCBIWWW.qblast(
    "blastn",  # program
    "nt",      # database
    record.seq
)

# Save
with open("result.xml", "w") as f:
    f.write(result_handle.read())

# Parse
with open("result.xml") as f:
    blast_records = NCBIXML.parse(f)
    
    for record in blast_records:
        for alignment in record.alignments:
            print(f"Hit: {alignment.title}")
            print(f"Length: {alignment.length}")
            
            for hsp in alignment.hsps:
                print(f"  E-value: {hsp.expect}")
                print(f"  Score: {hsp.score}")
                print(f"  Identity: {hsp.identities}/{hsp.align_length}")
                print(f"  Query: {hsp.query}")
                print(f"  Match: {hsp.match}")
                print(f"  Sbjct: {hsp.sbjct}")
""")

print("\n=== Local BLAST ===")
print("""
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio.Blast import NCBIXML

# Run local BLAST
blastn_cline = NcbiblastnCommandline(
    query="query.fasta",
    db="nt",
    evalue=1e-5,
    outfmt=5,  # XML
    out="result.xml"
)

stdout, stderr = blastn_cline()
print(stdout)
print(stderr)

# Parse XML
with open("result.xml") as f:
    blast_records = NCBIXML.parse(f)
    for record in blast_records:
        for alignment in record.alignments:
            for hsp in alignment.hsps:
                print(f"E-value: {hsp.expect}")
""")

print("\n=== Parse BLAST Output ===")
print("""
from Bio.Blast import NCBIXML

def parse_blast(xml_file):
    \"\"\"Parse BLAST XML output.\"\"\"
    with open(xml_file) as f:
        records = NCBIXML.parse(f)
        
        for record in records:
            print(f"Query: {record.query}")
            print(f"Length: {record.query_length}")
            
            for alignment in record.alignments:
                print(f"\\nHit: {alignment.title}")
                print(f"Length: {alignment.length}")
                
                for hsp in alignment.hsps:
                    print(f"  E-value: {hsp.expect}")
                    print(f"  Bit score: {hsp.score}")
                    print(f"  Identities: {hsp.identities}/{hsp.align_length} "
                          f"({100*hsp.identities/hsp.align_length:.1f}%)")
                    print(f"  Gaps: {hsp.gaps}")
                    print(f"  Strand: {hsp.strand}")
                    print(f"  Query: {hsp.query[:50]}...")
                    print(f"  Match: {hsp.match[:50]}...")
                    print(f"  Sbjct: {hsp.sbjct[:50]}...")
""")

print("\n=== Using Bio.SearchIO ===")
print("""
from Bio import SearchIO

# Parse (newer API)
for qresult in SearchIO.parse("result.xml", "blast-xml"):
    print(f"Query: {qresult.id}")
    print(f"Hits: {len(qresult)}")
    
    for hit in qresult:
        print(f"\\nHit: {hit.id}")
        print(f"Description: {hit.description}")
        print(f"E-value: {hit.evalue}")
        print(f"Bit score: {hit.bitscore}")
        
        for hsp in hit:
            print(f"  HSP: {hsp.query_start}-{hsp.query_end}")
            print(f"  Identity: {hsp.ident_num}/{hsp.aln_span}")
""")

print("\n=== Filter Results ===")
print("""
def filter_blast_results(xml_file, evalue_threshold=1e-5, 
                          min_identity=80, min_coverage=0.7):
    \"\"\"Filter BLAST results.\"\"\"
    filtered = []
    
    with open(xml_file) as f:
        records = NCBIXML.parse(f)
        
        for record in records:
            query_len = record.query_length
            
            for alignment in record.alignments:
                for hsp in alignment.hsps:
                    # Filters
                    if hsp.expect > evalue_threshold:
                        continue
                    
                    identity = hsp.identities / hsp.align_length
                    if identity < min_identity / 100:
                        continue
                    
                    coverage = hsp.align_length / query_len
                    if coverage < min_coverage:
                        continue
                    
                    filtered.append({
                        'hit': alignment.title,
                        'evalue': hsp.expect,
                        'bitscore': hsp.score,
                        'identity': identity,
                        'coverage': coverage,
                        'query': hsp.query,
                        'sbjct': hsp.sbjct,
                    })
    
    return filtered

# Use
results = filter_blast_results("result.xml", evalue_threshold=1e-10)
for r in results:
    print(f"{r['hit']}: E={r['evalue']}, Id={r['identity']:.2%}")
""")

print("\n=== Batch BLAST ===")
print("""
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio import SeqIO
import os

def batch_blast(input_fasta, output_dir, db="nt", evalue=1e-5):
    \"\"\"Run BLAST on multiple sequences.\"\"\"
    os.makedirs(output_dir, exist_ok=True)
    
    for record in SeqIO.parse(input_fasta, "fasta"):
        output = os.path.join(output_dir, f"{record.id}.xml")
        
        blastn = NcbiblastnCommandline(
            query="-",  # stdin
            db=db,
            evalue=evalue,
            outfmt=5,
            out=output
        )
        
        # Run
        stdout, stderr = blastn(stdin=str(record.seq))
        print(f"Processed: {record.id}")

# Use
# batch_blast("queries.fasta", "blast_results/", db="nt")
""")

print("\n=== Complete Example ===")
print("""
from Bio.Blast import NCBIWWW, NCBIXML
from Bio import SeqIO

def run_blast_remote(sequence, program="blastn", database="nt"):
    \"\"\"Run BLAST remotely and return parsed results.\"\"\"
    # Run
    result_handle = NCBIWWW.qblast(program, database, sequence)
    
    # Parse
    blast_records = NCBIXML.read(result_handle)
    
    # Extract
    results = []
    for alignment in blast_records.alignments:
        for hsp in alignment.hsps:
            results.append({
                'title': alignment.title,
                'length': alignment.length,
                'evalue': hsp.expect,
                'score': hsp.score,
                'identities': hsp.identities,
                'align_length': hsp.align_length,
                'query': hsp.query,
                'sbjct': hsp.sbjct,
            })
    
    return results

# Use
# seq = "ATGCGTACGTTAG" * 10
# results = run_blast_remote(seq)
# for r in results[:5]:
#     print(f"{r['title'][:60]} E={r['evalue']}")
""")

⌨️ Command-Line BLAST

🎯 BLAST from the Terminal

BLAST is typically run from the command line for local databases. The commands are flexible but have many options.

CLI

blastn: Nucleotide. blastp: Protein. Options: Query, db, evalue. Output: Format. Makeblastdb: Create DB. Batch: Loop.

⌨️ blastn⌨️ blastp⌨️ makeblastdb⌨️ outfmt

# Command-Line BLAST
print("=== Installation ===")
print("""
# Ubuntu/Debian
sudo apt-get install ncbi-blast+

# macOS
brew install blast

# Download from NCBI
# https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

# Check
blastn -version
""")

print("\n=== Create Database ===")
print("""
# Nucleotide database
makeblastdb -in sequences.fasta -dbtype nucl -out mydb

# Protein database
makeblastdb -in proteins.fasta -dbtype prot -out mydb

# Options
-in      Input FASTA
-dbtype  nucl or prot
-out     Database name
-title   Description
""")

print("\n=== Basic BLAST ===")
print("""
# blastn (DNA vs DNA)
blastn -query query.fasta -db nt -out results.txt

# blastp (protein vs protein)
blastp -query protein.fasta -db nr -out results.txt

# blastx (DNA vs protein)
blastx -query dna.fasta -db nr -out results.txt

# tblastn (protein vs DNA)
tblastn -query protein.fasta -db nt -out results.txt

# tblastx (DNA vs DNA, translated)
tblastx -query dna.fasta -db nt -out results.txt
""")

print("\n=== Common Options ===")
print("""
-query        Query file
-db           Database
-out          Output file
-evalue       E-value threshold (default 10)
-outfmt       Output format
-num_threads  CPU threads
-max_target_seqs  Max hits (default 500)
-word_size    Word size (11 for blastn, 3 for blastp)
-gapopen      Gap open penalty
-gapextend    Gap extend penalty
-matrix       Scoring matrix (BLOSUM62)
-strand       plus, minus, or both
-query_gencode  Genetic code
""")

print("\n=== Output Format ===")
print("""
# Pairwise (default)
blastn -query q.fasta -db nt -out results.txt

# Tabular
blastn -query q.fasta -db nt -outfmt 6 -out results.tsv

# Custom tabular
blastn -query q.fasta -db nt \\
    -outfmt "6 qseqid sseqid pident length evalue bitscore" \\
    -out results.tsv

# XML
blastn -query q.fasta -db nt -outfmt 5 -out results.xml

# JSON
blastn -query q.fasta -db nt -outfmt 15 -out results.json
""")

print("\n=== Example Commands ===")
print("""
# Basic blastn
blastn -query gene.fasta -db nt -out result.txt

# With E-value
blastn -query gene.fasta -db nt -evalue 1e-10 -out result.txt

# Multi-threaded
blastn -query gene.fasta -db nt -num_threads 8 -out result.txt

# Tabular output
blastn -query gene.fasta -db nt -outfmt 6 -out result.tsv

# Custom tabular
blastn -query gene.fasta -db nt \\
    -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \\
    -out result.tsv

# Limited hits
blastn -query gene.fasta -db nt -max_target_seqs 10 -out result.txt

# Both strands
blastn -query gene.fasta -db nt -strand both -out result.txt

# Blastp with BLOSUM62
blastp -query protein.fasta -db nr -matrix BLOSUM62 -out result.txt

# Custom gap penalties
blastp -query protein.fasta -db nr -gapopen 11 -gapextend 1 -out result.txt
""")

print("\n=== Batch Processing ===")
print("""
#!/bin/bash

# Loop over queries
for query in queries/*.fasta; do
    name=$(basename $query .fasta)
    blastn -query $query -db nt -out results/${name}.txt
done

# Parallel
ls queries/*.fasta | parallel -j 8 \\
    'blastn -query {} -db nt -out results/{/.}.txt'
""")

print("\n=== Making Database ===")
print("""
# Download sequences
wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.*.tar.gz
tar -xzf nt.*.tar.gz

# Or make custom
makeblastdb -in my_sequences.fasta -dbtype nucl -out mydb

# Check
blastdbcmd -db mydb -info
blastdbcmd -db mydb -entry all -outfmt %f -out all.fasta
""")

print("\n=== Tips ===")
print("""
1. Use -num_threads for speed
2. Precompute database
3. Use tabular output for parsing
4. Set evalue threshold appropriately
5. Consider database size
6. Use -max_target_seqs to limit
7. Check strand for DNA
8. Save command for reproducibility
""")

⚖️ vs Smith-Waterman

🎯 BLAST vs Smith-Waterman

BLAST is a heuristic (fast, may miss); Smith-Waterman is exact (slow, optimal). They serve different purposes.

Comparison

BLAST: Heuristic, fast. SW: Exact, slow. BLAST: 100-1000x faster. SW: Guaranteed optimal. BLAST: Database search. SW: Confirmation.

⚖️ BLAST⚖️ Smith-Waterman⚖️ Heuristic⚖️ Exact

Aspect	BLAST	Smith-Waterman
Type	Heuristic	Exact
Optimality	Not guaranteed	Guaranteed optimal
Speed	Fast (seconds)	Slow (minutes-hours)
Complexity	O(m + n × k)	O(mn)
Sensitivity	High (but may miss)	Highest
Speed vs SW	100-1000x faster	1x
Algorithm	Seed-and-extend	Dynamic programming
Statistics	E-value, bit score	Score (raw)
Use case	Database search	Confirmation, sensitive
When to use	Millions of queries	Few, high-value queries
Tools	BLAST+, DIAMOND	SSEARCH, SSEARCH, parasail

# When to Use Which
print("=== Use BLAST WHEN ===")
print("""
✓ Database search (millions of sequences)
✓ Speed is critical
✓ Initial screening
✓ Most cases in practice
✓ Millions of queries

Examples:
  - Annotating a genome
  - Finding homologs
  - Metagenomics
  - Routine analysis
""")

print("\n=== Use Smith-Waterman WHEN ===")
print("""
✓ Confirm a BLAST hit
✓ Sensitive search for remote homologs
✓ Small number of queries
✓ Optimal alignment needed
✓ BLAST missed something

Examples:
  - Confirming distant homology
  - Validating a hit
  - High-value targets
  - Specialized analysis
""")

print("\n=== Combined Workflow ===")
print("""
Two-step approach:
  1. BLAST (fast, heuristic)
     - Find candidates
     - Filter by E-value
  2. Smith-Waterman (slow, exact)
     - Confirm candidates
     - Get optimal alignment
     - Verify significance

This is standard practice:
  - Fast screening
  - Accurate confirmation
  - Best of both
""")

print("\n=== Speed Comparison ===")
print(f"{'Method':<20} {'Time (1M seqs)':<20} {'Sensitivity':<15}")
print("-" * 55)
print(f"{'BLAST':<20} {'~1 minute':<20} {'High':<15}")
print(f"{'DIAMOND':<20} {'~10 seconds':<20} {'Medium':<15}")
print(f"{'Smith-Waterman':<20} {'~10 hours':<20} {'Highest':<15}")
print(f"{'GPU SW':<20} {'~6 minutes':<20} {'Highest':<15}")

💼 Use Cases

🎯 Where BLAST Is Used

BLAST is used in virtually every area of molecular biology: gene annotation, metagenomics, phylogenetics, drug discovery, and more.

Use Cases

Genome annotation: Identify genes. Homology: Find relatives. Metagenomics: Classify species. Phylogenetics: Evolutionary distance. Drug: Target identification. Diagnostics: Pathogen detection.

💼 Genome💼 Homology💼 Metagenomics💼 Phylogenetics

# Use Cases
print("=== 1. Genome Annotation ===")
print("""
Use case:
  - Identify genes in new genome
  - Predict function
  - Annotate ORFs

Workflow:
  1. Sequence genome
  2. Predict ORFs
  3. BLAST each ORF vs nr
  4. Assign function based on hits

Tools:
  - Prokka (prokaryotes)
  - MAKER (eukaryotes)
  - NCBI PGAAP
""")

print("\n=== 2. Homology Detection ===")
print("""
Use case:
  - Find homologs
  - Determine relationships
  - Identify conserved regions

Workflow:
  1. Query sequence
  2. BLAST vs database
  3. Filter by E-value
  4. Analyze top hits
  5. Build phylogeny

Example:
  - Novel protein
  - BLAST finds similar proteins
  - Function inferred
""")

print("\n=== 3. Metagenomics ===")
print("""
Use case:
  - Classify environmental samples
  - Identify species
  - Functional profiling

Workflow:
  1. Sequence sample
  2. BLAST reads vs reference
  3. Classify by top hit
  4. Quantify abundance

Tools:
  - Kraken (k-mer based)
  - MetaPhlAn (marker genes)
  - DIAMOND (fast BLAST)

Example:
  - Gut microbiome
  - Ocean samples
  - Soil samples
""")

print("\n=== 4. Phylogenetics ===")
print("""
Use case:
  - Evolutionary relationships
  - Tree construction
  - Species comparison

Workflow:
  1. Collect sequences
  2. BLAST to find homologs
  3. Multiple alignment
  4. Build tree

Example:
  - 16S rRNA
  - Cytochrome c
  - Mitochondrial DNA

Tools:
  - MEGA
  - PhyML
  - RAxML
""")

print("\n=== 5. Variant Detection ===")
print("""
Use case:
  - Find SNPs
  - Detect disease variants
  - Compare to reference

Workflow:
  1. Sequence patient DNA
  2. BLAST vs reference
  3. Identify variants
  4. Annotate impact

Example:
  - Cancer genomics
  - Rare disease
  - Pharmacogenomics
""")

print("\n=== 6. Drug Discovery ===")
print("""
Use case:
  - Target identification
  - Off-target prediction
  - Drug repurposing

Workflow:
  1. Target protein
  2. BLAST vs proteome
  3. Find similar proteins
  4. Predict off-targets

Example:
  - Kinase inhibitors
  - GPCR drugs
  - Antibiotics
""")

print("\n=== 7. Diagnostics ===")
print("""
Use case:
  - Pathogen detection
  - Outbreak tracking
  - AMR prediction

Workflow:
  1. Sequence sample
  2. BLAST vs pathogen DB
  3. Identify species
  4. Detect resistance

Example:
  - COVID-19 variants
  - Tuberculosis
  - Antimicrobial resistance
""")

print("\n=== 8. CRISPR Design ===")
print("""
Use case:
  - Off-target prediction
  - Guide RNA design
  - Specificity check

Workflow:
  1. Candidate guide
  2. BLAST vs genome
  3. Find off-targets
  4. Rank guides

Tools:
  - Cas-OFFinder
  - Benchling
  - CRISPRscan
""")

print("\n=== Use Case Comparison ===")
print(f"{'Use Case':<25} {'Variant':<15} {'Database':<20}")
print("-" * 60)
print(f"{'Genome annotation':<25} {'blastx':<15} {'nr':<20}")
print(f"{'Homology':<25} {'blastp':<15} {'nr':<20}")
print(f"{'Metagenomics':<25} {'blastx':<15} {'nt':<20}")
print(f"{'Phylogenetics':<25} {'blastn':<15} {'nt':<20}")
print(f"{'Variant detection':<25} {'blastn':<15} {'ref':<20}")
print(f"{'Drug discovery':<25} {'blastp':<15} {'pdb':<20}")
print(f"{'Diagnostics':<25} {'blastn':<15} {'pathogen':<20}")
print(f"{'CRISPR':<25} {'blastn':<15} {'genome':<20}")

❓ FAQ

🎯 Frequently Asked Questions

Answers to common questions about BLAST.

What is BLAST?

BLAST (Basic Local Alignment Search Tool) is a heuristic algorithm for comparing biological sequences (DNA, RNA, or protein) against sequence databases. It finds regions of local similarity by using a seed-and-extend strategy, trading some sensitivity for dramatic speed improvements over exact algorithms like Smith-Waterman.

What does E-value mean in BLAST?

The E-value (Expectation value) is the number of hits with a score at least as good as the observed score that would be expected by chance when searching a database of a particular size. Smaller E-values indicate more significant matches. E-values below 1e-5 are typically considered significant; below 1e-50 are highly significant.

What is the difference between BLAST and Smith-Waterman?

BLAST is a heuristic algorithm that uses seed-and-extend to find candidate alignments quickly — it's 100-1000x faster than Smith-Waterman. Smith-Waterman is exact and guaranteed optimal but much slower. BLAST may miss some remote homologs that Smith-Waterman would find. Use BLAST for database searches; use Smith-Waterman for confirmation and high-sensitivity alignment.

What are the 5 BLAST variants?

The five main BLAST programs are: (1) blastn — nucleotide query vs nucleotide database; (2) blastp — protein query vs protein database; (3) blastx — translated nucleotide query vs protein database; (4) tblastn — protein query vs translated nucleotide database; (5) tblastx — translated nucleotide query vs translated nucleotide database.

What is a bit score in BLAST?

The bit score (S') is a normalized alignment score that is independent of database size. It is derived from the raw alignment score (S) using the statistical parameters of the scoring system: S' = (λS - ln K) / ln 2. Higher bit scores indicate better alignments. Bit scores can be compared across different database searches, unlike raw scores.

What is the best E-value threshold?

The default is 10 (very permissive). For most analyses, use 1e-5 or stricter. For high-confidence homologs, use 1e-20. The choice depends on your tolerance for false positives vs false negatives. For database annotation, stricter thresholds (1e-10 or lower) are common.

Can BLAST find distantly related sequences?

Sometimes. BLAST is less sensitive than Smith-Waterman for distant homologs. To improve sensitivity: use smaller word size, lower E-value threshold, BLOSUM45 matrix, and consider PSI-BLAST (iterative). For very distant homologs, use HMMER (profile HMM), MMseqs2, or DIAMOND in sensitive mode.

How do I run BLAST locally?

Install BLAST+ from NCBI (apt, brew, or download). Download a database (e.g., nt, nr) or create your own with makeblastdb. Then run: blastn -query query.fasta -db nt -out results.txt. Use -num_threads for speed and -outfmt 6 for tabular output.

What is the difference between blastn and blastp?

blastn compares nucleotide sequences directly (DNA vs DNA). blastp compares protein sequences directly (protein vs protein). blastp is usually more sensitive for detecting homology because proteins are more conserved than DNA (due to the genetic code's redundancy).

What are common BLAST alternatives?

DIAMOND (100-1000x faster, protein), MMseqs2 (very fast, sensitive), PSI-BLAST (iterative, sensitive), HMMER (profile HMM, most sensitive for protein families), SSEARCH (Smith-Waterman, exact), BLASTP+ (NCBI's newer version).

# Quick Reference
print("=== Quick Answers ===")
print("""
Q: What is BLAST?
A: Heuristic sequence database search.

Q: What does E-value mean?
A: Expected random hits; smaller = better.

Q: Best threshold?
A: 1e-5 typical, 1e-10 strict.

Q: BLAST vs SW?
A: BLAST heuristic fast, SW exact slow.

Q: Best variant for genes?
A: blastx (DNA vs protein DB).

Q: Best variant for proteins?
A: blastp.

Q: How to speed up?
A: Use DIAMOND or MMseqs2.

Q: How to increase sensitivity?
A: PSI-BLAST or HMMER.

Q: Best database?
A: nr (protein), nt (nucleotide).

Q: Free?
A: Yes — NCBI BLAST+.
""")

🚀 Future

🎯 Where BLAST Is Going

The future is AI-augmented, hardware-accelerated, cloud-native, and integrated with pangenomics and long reads.

Future Directions

AI: Learned scoring. Hardware: GPU, FPGA. Cloud: Serverless. Pangenome: Graph. Long reads: Adaptive. Integration: Unified.

🚀 AI🚀 Hardware🚀 Cloud🚀 Pangenome

# Future of BLAST
print("=== 1. AI-Augmented ===")
print("""
Trend: ML for alignment

Ideas:
  - Learned scoring matrices
  - Neural alignment
  - Transformer-based
  - Protein language models

Examples:
  - ESM (Meta)
  - AlphaFold (DeepMind)
  - DNABERT
  - Nucleotide Transformer

Future:
  - Learned substitution
  - Better sensitivity
  - Faster
""")

print("\n=== 2. Hardware Acceleration ===")
print("""
Trend: Faster BLAST

Current:
  - CPU
  - SIMD (SSE, AVX)
  - 10x speedup

Future:
  - GPU
  - FPGA
  - ASIC
  - 100-1000x speedup

Examples:
  - NVIDIA Parabricks
  - Illumina DRAGEN
  - AWS F1 (FPGA)
  - Cerebras
""")

print("\n=== 3. Cloud-Native ===")
print("""
Trend: Cloud BLAST

Current:
  - Self-managed servers
  - Manual scaling

Future:
  - Serverless BLAST
  - Auto-scaling
  - Pay per use
  - Global

Examples:
  - AWS HealthOmics
  - Azure Genomics
  - Google Cloud Life Sciences
  - Terra
  - Galaxy
""")

print("\n=== 4. Pangenomics ===")
print("""
Trend: Graph genomes

Traditional:
  - Linear reference
  - BLAST aligns to one sequence

Pangenome:
  - Graph reference
  - Multiple genomes
  - Better for diversity

Algorithms:
  - Graph BLAST
  - Variation graphs
  - Extension of seed-and-extend

Tools:
  - vg toolkit
  - GraphAligner
  - Minigraph
""")

print("\n=== 5. Long Reads ===")
print("""
Trend: Long-read sequencing

Current:
  - Short reads (100-300 bp)
  - Standard BLAST

Future:
  - Long reads (10-100 kb)
  - Adaptive alignment
  - Banded BLAST
  - Sparse seed

Tools:
  - minimap2
  - GraphMap
  - deSALT
  - BLAST with long-read mode
""")

print("\n=== 6. Integration ===")
print("""
Trend: Unified platforms

Current:
  - BLAST separate
  - Many tools

Future:
  - Integrated pipelines
  - Workflow managers
  - Cloud platforms
  - Reproducible

Examples:
  - Galaxy
  - Nextflow
  - Snakemake
  - Terra
""")

# TIMELINE
print("\n=== Timeline ===")
print("""
Past:
  1990: BLAST
  1997: PSI-BLAST
  2000: BLAST+ (NCBI)
  2015: DIAMOND
  2017: MMseqs2
  2020: GPU BLAST

Present:
  2025: Hardware-accelerated
  2026: AI-augmented

Future:
  2030: Serverless?
  2035: Pangenome-native?
  2040: Quantum?
""")

print("\n=== Final Thoughts ===")
print("""
BLAST is:
  - Heuristic sequence search
  - Fast and scalable
  - Ubiquitous in biology
  - Foundation of genomics
  - Still essential

The future is:
  - AI-augmented
  - Hardware-accelerated
  - Cloud-native
  - Pangenome
  - Long reads

"BLAST is the Google of biology."

Keep exploring. Keep building.
""")