Searching databases with BLAST

Seeding, extending, and the price of speed

From last time

Smith-Waterman finds the best local alignment: The best match between part of one sequence and part of another.

TTTTACGTCA
    ||||||
    ACGTCAGGGG

“Is some region of this database sequence similar to some region of my sequence?” is a local alignment question.

That is the question you ask every time you run BLAST.

When you use BLAST

  • Paste in a sequence and click a button
  • A few seconds later: A ranked list of similar sequences
  • Searched against GenBank: About 60 trillion bases1

How is that possible?!

Can we just use Smith-Waterman?

The scoring table has \(m \times n\) cells (\(O(mn)\)).

One 150-base read against GenBank:

\[ 150 \times (6 \times 10^{13}) \approx 9 \times 10^{15} \text{ cells} \]

  • Optimized code, a billion cells per second: Over 3 months
  • Our Python code: A century or more (you will measure this today)
  • One sequencing run has billions of reads
    • Click, go on summer vacation, repeat

A faster computer does not fix this. We need to avoid doing most of the work.

Heuristic algorithms

From our algorithms lecture:

  • Smith-Waterman is optimal but far too slow to scan all of GenBank
  • BLAST first finds short exact seed matches, then only extends around those

“Approximate answer today” beats “exact answer in nine years.”

Today we unpack that second bullet: BLAST1, the Basic Local Alignment Search Tool.

Seed and extend

The key insight

Two homologous sequences almost always share at least one short, exact match.

TGCTTACCAGTCAGGT
TGCTTACGAGTCATTG

So:

  1. Find short exact matches (seeds) — fast!
  2. Do alignment work only near the seeds
  3. Ignore the rest of the database

Step 0: Index the database, once

Finding exact matches of length \(k\) fast?

  • We have done this!

The \(k\)-mer dictionary from our string search lecture:

{"TGCT": [7], "GCTT": [8], "CTTA": [9], ...}
  • Build it once for the whole database
  • Every lookup after that is \(O(1)\)
  • NCBI builds its BLAST databases long before you submit a query

Step 1: Seeds

Look up every \(k\)-mer of the query (\(k = 4\) below):

query:   TGCTTACCAGTCAGGT
subject: GGATCCATGCTTACGAGTCATTGCAACGT
Query \(k\)-mer Query position Subject position(s)
TGCT 0 7
GCTT 1 8
CTTA 2 9
TTAC 3 10
TACC 4 —
…
AGTC 8 15
GTCA 9 16

Homologous seeds often line up on diagonals

Seed Query pos Subject pos Subject \(-\) query
TGCT 0 7 7
GCTT 1 8 7
CTTA 2 9 7
TTAC 3 10 7
AGTC 8 15 7
GTCA 9 16 7

All six are on the same diagonal of a dotplot: Pieces of the same alignment!

A seed alone on its own diagonal is usually just chance.

Step 2: Extend without gaps

Grow the seed one base at a time, keeping a running score (match \(+1\), mismatch \(-1\)):

seed
query TGCT T A C C A G T C A G G T
subject TGCT T A C G A G T C A T T G
score 4 5 6 7 6 7 8 9 10 11 10 9 8

X-drop rule: stop once the score falls more than \(X\) below the best score seen (here \(X = 2\)), then trim back to the best.

The result is a high-scoring segment pair (HSP) with a score of 11.

Indels break ungapped extension

Alpha: CCATGCTAT------CTCTGGGA
Wuhan: CCATGCTATACATGTCTCTGGGA
  • Seeds on both sides of the indel, but on two different diagonals
  • Ungapped extension gives two separate HSPs

Gapped BLAST1: Run a gapped alignment, but only in the short region around a good HSP.

Run the Smith-Waterman aligner on ~200 bases, NOT 60 trillion.

Step 3: Is the hit meaningful?

In a database of 60 trillion bases, many matches happen by chance.

BLAST reports an E-value1: The number of hits scoring at least this well that you would expect by chance in a database this big.

\[ E = K\,m\,n\,e^{-\lambda S} \]

  • Bigger query (\(m\)) or database (\(n\)): More chance hits
  • Higher score (\(S\)): Exponentially fewer chance hits
  • \(E = 10^{-50}\): Essentially impossible by chance; \(E = 5\): Probably just chance

The price of speed

Missing homologous sequences

BLAST only looks where it finds a seed.

A homolog with no exact match of \(k\) consecutive bases is invisible to BLAST, no matter how good the full alignment would be.

TGCTTACCAGTCAGGT
| || ||| || || |
TCCTAACCTGTAAGCT

11 of 16 bases identical, but never more than 3 in a row

  • With \(k = 4\), we miss it!

The word size (\(k\)) trade-off

Number of chance seeds for a 150-base query against GenBank (\(\approx m\,n / 4^k\)):

Word size \(k\) Chance seeds Needs an exact run of
7 \(\sim 5 \times 10^{11}\) 7 bases
11 \(\sim 2 \times 10^{9}\) 11 bases
28 \(\sim 0.1\) 28 bases

Small \(k\)

More sensitive, much slower

Large \(k\)

Much faster, misses divergent homologs

NCBI’s tool choices

Program (task) Word size Designed for
megablast (default) 28 Highly similar sequences
blastn 11 Somewhat similar sequences
blastn-short 7 Short queries (e.g., primers)
Protein BLAST 3 Proteins (20-letter alphabet)

When you pick an option on the BLAST web page, you are choosing a point along the speed/sensitivity trade-off continuum!

Applying what we learned

Which BLAST program would you use?

  • Which species did this 16S sequence from a soil sample come from?
  • Does a human genome contain a homolog of this frog gene?
  • Where does this 20-base PCR primer bind in the human genome?

BLAST does a bit more

  • Searches both strands of DNA
  • Masks low-complexity sequences (e.g., ATATATATAT) that would seed everywhere
  • Protein BLAST requires two seeds near each other before extending
  • Reports multiple HSPs per database sequence

But the core is what you will build today: Seed, extend, evaluate

Summary

  • Optimal local alignment is \(O(mn)\), but \(n\) is tens of trillions!
  • Seed: Find short exact matches with an index built once
  • Extend: Ungapped with the X-drop rule, then gapped only near good HSPs
  • Evaluate: E-values suggest whether a hit is better than chance
  • Word size sets the speed/sensitivity trade-off
  • BLAST gives up the guarantee of optimality to get an answer fast

In the exercise

You will build BLAST from pieces you mostly already have:

  1. Turn your Needleman-Wunsch code into Smith-Waterman (three small changes)
  2. Time an optimal search of five SARS-CoV-2 genomes
  3. find_seeds: Find seeds with your \(k\)-mer dictionary
  4. Use seed-and-extend to identify a mystery read
  5. Explore the word size trade-off

Watch for

Which variant the mystery read came from, and which variant is its runner-up. The sequence alignment exercise holds the clue.

Exercise

  1. Log into molab
  2. Open the exercise in molab
  3. Click “Fork” button

References

Altschul, Stephen F., Warren Gish, Webb Miller, Eugene W. Myers, and David J. Lipman. 1990. “Basic Local Alignment Search Tool.” Journal of Molecular Biology 215 (3): 403–10. https://doi.org/10.1016/S0022-2836(05)80360-2.
Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schäffer, et al. 1997. “Gapped BLAST and PSI-BLAST: A New Generation of Protein Database Search Programs.” Nucleic Acids Research 25 (17): 3389–402. https://doi.org/10.1093/nar/25.17.3389.
Karlin, Samuel, and Stephen F. Altschul. 1990. “Methods for Assessing the Statistical Significance of Molecular Sequence Features by Using General Scoring Schemes.” Proceedings of the National Academy of Sciences 87 (6): 2264–68. https://doi.org/10.1073/pnas.87.6.2264.

Acknowledgements

I used Anthropic’s Claude (Opus 5 model) to generate a draft of these slides.