Generated by the AI hub at UTK, beware!

Module goal: Understand how we search enormous sequence databases quickly, and why "quickly" required giving up on guaranteed-optimal alignment.

You have a sequence. Somewhere in a database of billions of sequences, there may be something related to it. How do you find it before you retire?

The mathematically correct answer is dynamic programming — Smith-Waterman for local alignment. It is guaranteed to find the optimal alignment between two sequences. It is also $$O(mn)$$: the work scales with the length of your query times the length of the database.

$$\text{time} \propto (\text{query length}) \times (\text{database length})$$

For a 300-residue protein against a 200-billion-residue database, that's $$6 \times 10^{13}$$ cell computations per query. Correct, and useless.

The central trade-off of this entire module: Every tool we discuss buys speed by giving up the guarantee of optimality. The interesting question is never "is it exact?" (it isn't) but "what does it miss, and do I care?"

Basic Local Alignment Search Tool — the heuristic that made database search practical.

Real homologs almost always share at least one short stretch of near-identical sequence. So don't align everything. Find those short stretches first, and only do expensive work around them.

This is called seed-and-extend.

BLAST chops your query into overlapping short substrings called words (default length $$w = 3$$ for proteins, 11 for nucleotides).

QUERY:   M K V L S T A E G
          \_/
          MKV
           \_/
           KVL
            \_/
            VLS      ... and so on

Critically, for proteins BLAST doesn't only look for exact word matches. It builds a neighborhood of every word that scores above a threshold $$T$$ against the query word, using a substitution matrix (BLOSUM62 by default).

Query word:  MKV

Neighborhood words (score ≥ T under BLOSUM62):
   MKV  (exact)
   MRV  (K→R is a conservative swap)
   LKV  (M→L is conservative)
   IKV
   ...

The database is pre-indexed, so BLAST can ask "where does any of these words appear?" and get an answer instantly. Each location is a seed — a candidate worth investigating.

Why the neighborhood matters: proteins diverge. Requiring exact 3-mer matches would miss most remote homologs. Requiring exact 11-mers in DNA is fine because DNA homology decays differently. This is why protein search is more sensitive than DNA search for distant relationships — and why DIAMOND (Section 3) works on protein space.

Once BLAST has a seed, it extends the match outward in both directions, one residue at a time, accumulating a score.

                  seed
                   ↓
  query   ... A E L M K V L S T A E ...
              | |   | | | | |   |
  subject ... A D L M K V L S G A E ...
              ←──── extend ────→

  score:  ..12 14 13 18 22 26 31 30 24 ..
                              ↑
                          MAXIMUM — stop and trim back here

Extension stops when the score drops more than $$X$$ below the best score seen so far (the X-drop rule), then trims back to the maximum. The surviving segment is a High-Scoring Segment Pair (HSP).

An HSP is: a pair of equal-length subsequences (one from query, one from subject) whose alignment score cannot be improved by extending or shortening it.

⚠️ Note the 1990 limitation: HSPs are ungapped. No insertions, no deletions. A real homolog with an indel in the middle shows up as two separate HSPs rather than one alignment. Fixing this is what 1997 is about.

This is the piece that made BLAST scientifically useful rather than just fast. Anyone can report "here are some matches." BLAST tells you which matches would have happened by chance anyway.

The E-value (Expect value) answers:

In a database this size, how many alignments this good or better would I expect to find purely by chance?

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

How to read it:

1. An E-value is not a p-value. It's an expected count, not a probability. It can exceed 1. (They converge for small values: $$p = 1 - e^{-E} \approx E$$ when $$E \ll 1$$.)

2. E-values depend on database size. The same alignment gets a worse E-value in a bigger database, because more chances to find it by luck. Your hit from three years ago may not clear the same threshold today.

$$E \propto n \quad \Rightarrow \quad \text{database doubles} \Rightarrow \text{E-value doubles}$$

3. Low E-value ≠ same function. It means unlikely to be chance similarity. Homology is a statement about shared ancestry, not shared job description.

Optimal local alignment scores between random sequences follow a Gumbel (extreme value) distribution, not a normal one — because you're taking a maximum over many possible alignments, and maxima of many draws have their own characteristic distribution with a long right tail. That heavy tail is exactly why an eyeball assessment of "that looks like a good score" is unreliable and the formal statistics matter.

   frequency
      │      ╭─╮
      │     ╱   ╲
      │    ╱     ╲___
      │   ╱          ╲______        ← long right tail:
      │  ╱                   ╲____    high scores DO occur
      │ ╱                         ╲__  by chance
      └────────────────────────────────  alignment score

Two upgrades, both still in daily use.

The two-hit requirement. Instead of extending from every single word hit, BLAST 2.0 only extends when it finds two non-overlapping word hits on the same diagonal within distance $$A$$ of each other.

DIAGONAL VIEW (query vs. subject offsets)

  ─────●──────────────●─────────────  ← two hits, close together: EXTEND
        ← distance A →

  ─────●───────────────────────────●─  ← too far apart: skip
  
  ──────────●──────────────────────── ← lone hit: skip

Counterintuitively, this makes BLAST both faster (far fewer extensions attempted) and lets them lower the threshold $$T$$ (more seeds generated, better sensitivity). Fewer, better-targeted extensions.

Gapped extension. Once a promising ungapped HSP is found, BLAST runs a restricted dynamic-programming alignment around it — allowing insertions and deletions, but only exploring a narrow band of the DP matrix instead of the whole thing.

Ungapped (1990):          Gapped (1997):

QUERY  AELMKVL            QUERY  AELMKVLSTAEG
       |||||||                   |||||||   ||||
SUBJ   AELMKVL            SUBJ   AELMKVL---AEG
       (one HSP)                 (one alignment, gap included)

...and separately         Biologically this is ONE homolog
QUERY  AEG                with a 3-residue deletion, and now
       |||                it's scored as one.
SUBJ   AEG
       (second HSP)

The practical result: one alignment per homolog instead of a scattered set of fragments, and a score that reflects the full relationship.

This is the conceptually bigger idea.

The problem: BLOSUM62 is a generic substitution matrix. It says "in proteins on average, K→R is a tolerable swap." But it knows nothing about your protein family. A catalytic cysteine that can never change is scored the same as a surface cysteine that changes freely.

A Position-Specific Scoring Matrix (PSSM) replaces the one-size-fits-all matrix with a per-position model of what's allowed at that position in this family.

GENERIC (BLOSUM62)           POSITION-SPECIFIC (PSSM)
one matrix, all positions    a column of scores per position

           A  C  D ...        pos1  pos2  pos3  pos4  pos5
        A  4  0 -2            A  -1    2   -3     5   -2
        C  0  9 -3            C  -4   -1   +9*   -4   -1
        D -2 -3  6            D   3   -2   -5    -1    6
        ...                   ...
                                        ↑
                              "at position 3, ONLY cysteine.
                               This is the catalytic residue."

PSI-BLAST (Position-Specific Iterated BLAST) builds this matrix automatically through iteration:

flowchart TD
    A["Round 1: standard BLAST<br/>with BLOSUM62"] --> B["Collect hits above<br/>threshold (e.g. E &lt; 0.005)"]
    B --> C["Build multiple alignment<br/>of those hits"]
    C --> D["Derive PSSM:<br/>which residues are<br/>conserved at each position?"]
    D --> E["Round 2+: search database<br/>again using the PSSM<br/>instead of BLOSUM62"]
    E --> F{"New hits found?"}
    F -->|Yes| B
    F -->|No| G["Converged.<br/>Report results."]
    
    style D fill:#e8f0fe,stroke:#3367d6
    style G fill:#e6f4ea,stroke:#137333

Why this is powerful: the profile progressively learns the family's signature and can detect homologs at sequence identities far below what standard BLAST reaches — the "twilight zone" around 20–30% identity where pairwise methods fail.

Why this is dangerous — profile drift. If one false positive enters the hit set in round 2, it contributes its residues to the PSSM in round 3, which then pulls in its relatives in round 4. The profile wanders away from your original query.

Practical rule: Use a strict inclusion threshold. Inspect what gets added each round. Don't run PSI-BLAST to convergence unattended and trust the output.

Same seed-and-extend philosophy, re-engineered for the sequencing era.

By 2015 the typical task wasn't "search one protein." It was "take 100 million short reads from a metagenome and assign each one a functional annotation." That means 100 million BLASTX searches (translating DNA reads in 6 frames and searching protein space).

BLASTX on that workload: weeks to months of CPU time.

How BLAST does it: index the database once. Then for each query, walk through the query's words and look each one up.

The hidden cost: each lookup jumps to an unpredictable place in a multi-gigabyte index. That's a cache miss — the CPU stalls waiting on RAM. When you do this billions of times, you spend most of your runtime waiting on memory, not computing.

How DIAMOND does it: index both the query set and the database. Then sort both lists of seeds and walk through them together, in order — like merging two sorted lists.

BLAST (single index):              DIAMOND (double index):

for each query word:               sorted query seeds:  sorted DB seeds:
   → random jump into                  AAC  ────────────── AAC
     the DB index                      AAC  ────────────── AAC
   → cache miss                        AAG  ────────────── AAG
   → CPU stalls                        ACD  ────────────── ACD
                                       ACF  ────────────── ACF
[unpredictable memory access]      
                                   [linear sweep, both lists
                                    read front to back —
                                    the CPU can prefetch]

Because both lists are sorted, memory access becomes sequential and predictable. Modern CPUs are enormously good at sequential access and terrible at random access. DIAMOND is not doing less algorithmic work — it's doing the same work in an order the hardware likes.

Additional design choices:

• Reduced amino acid alphabet — group biochemically similar residues so seeds tolerate conservative substitutions without exploding the neighborhood
• Spaced seeds — patterns like 11*1**1*1**11*111 that require matches only at specified positions, catching diverged homologs that contiguous seeds miss

Discussion hook: 20,000× is not a normal software speedup. It's the difference between a computation you plan a grant around and one you run over lunch. Speedups of that magnitude don't just accelerate existing science — they change which questions are askable.

DIAMOND 2015 was tuned for short reads. For full-length proteins and remote homology, it lost real sensitivity relative to BLASTP. If you needed BLASTP-quality answers, you still ran BLASTP.

New sensitivity modes let you dial the speed/sensitivity trade-off explicitly:

FAST ◄─────────────────────────────────────────────► SENSITIVE

--fast   (default)   --sensitive   --very-sensitive   --ultra-sensitive
   │                                      │                  │
   │                                      └──────────────────┘
   │                                       BLASTP-equivalent sensitivity
   │                                       at roughly 80× the speed
   │
   └─ tuned for short reads, maximum throughput

The headline: --very-sensitive and --ultra-sensitive match BLASTP sensitivity while running approximately 80× faster. The trade-off is now yours to make explicitly, rather than being forced on you by the tool.

The paper motivates this with the Earth BioGenome Project — the effort to sequence ~1.5 million eukaryotic species. Comparing all those proteomes against each other is an all-vs-all problem, and all-vs-all scales quadratically:

$$\text{comparisons} \propto N^2$$

At that scale, BLAST's runtime isn't inconvenient — the analysis is simply not performable. The paper estimates months of supercomputer time for tasks DIAMOND completes in hours.

Connect this back to week one. Stephens et al. argued genomics is a "four-headed beast": acquisition, storage, distribution, and analysis. Sequencing cost has collapsed — acquisition is largely solved. DIAMOND v2 is what solving the analysis head looks like in practice. Algorithmic efficiency isn't a convenience; at these scales it determines what science exists.

timeline
    title Same core idea, four scales of data
    1990 : BLAST
         : Seed and extend
         : Ungapped HSPs
         : E-value statistics
    1997 : Gapped BLAST + PSI-BLAST
         : Two-hit seeding
         : Gapped alignment
         : PSSMs and iteration
    2015 : DIAMOND
         : Double indexing
         : Cache-friendly design
         : ~20,000x on short reads
    2021 : DIAMOND v2
         : Tunable sensitivity modes
         : BLASTP-equivalent at ~80x
         : Tree-of-life scale

The through-line: the conceptual algorithm barely changed in 30 years. What changed is the size of the problem, and each paper is a response to a bottleneck that only became visible at the new scale.

Whatever you run, a reproducible methods section states:

• Tool and version (DIAMOND v0.9 and v2 give different results)
• Database and the date/version you downloaded it — E-values depend on database size
• Sensitivity mode / all non-default parameters
• E-value threshold, and any coverage or identity filters
• For PSI-BLAST: number of iterations and inclusion threshold

1. You BLAST a sequence in 2019 and get $$E = 0.001$$. You rerun the identical search in 2026 against the current database and get $$E = 0.02$$. The alignment is byte-for-byte identical. What changed, and which result should you report?

2. Your protein has a 30-residue insertion relative to its closest homolog. Sketch what BLAST 1990 reports versus BLAST 1997. Which is easier to interpret, and why does it matter for annotation?

3. DIAMOND doesn't use a fundamentally better algorithm than BLAST — the asymptotic complexity is comparable. Where does the 20,000× actually come from? What does that suggest about how you should think about performance in your own analyses?

4. You run PSI-BLAST for 5 iterations. Between round 2 and round 3 the hit count jumps from 40 to 900. Is this a discovery or a failure? How would you tell the difference?

5. Two sequences share 95% identity over 12 residues. Two others share 28% identity over 400 residues. Which is more likely to be a genuine homolog, and which statistic tells you?

6. Explain, without using the words "p-value" or "significant," what $$E = 10^{-8}$$ means to a colleague who has never run BLAST.

• Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ (1990). Basic local alignment search tool. J Mol Biol 215(3):403–410. https://doi.org/10.1016/S0022-2836(05)80360-2
• Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25(17):3389–3402. https://doi.org/10.1093/nar/25.17.3389
• Buchfink B, Xie C, Huson DH (2015). Fast and sensitive protein alignment using DIAMOND. Nature Methods 12:59–60. https://doi.org/10.1038/nmeth.3176
• Buchfink B, Reuter K, Drost HG (2021). Sensitive protein alignments at tree-of-life scale using DIAMOND. Nature Methods 18:366–368. https://doi.org/10.1038/s41592-021-01101-x

• NCBI BLAST Handbook — https://www.ncbi.nlm.nih.gov/books/NBK153387/
• DIAMOND documentation — https://github.com/bbuchfink/diamond/wiki