Sequence Alignment Algorithms: How BLAST Actually Searches

BLAST is not an alignment algorithm with a speed trick attached. It is a filter that finds candidate locations cheaply and only then runs alignment on them, and almost every parameter you can set controls the filter rather than the alignment.

The dynamic program BLAST is avoiding

Optimal local alignment of two sequences is solved exactly by the Smith-Waterman algorithm, published by Temple Smith and Michael Waterman in 1981 as a modification of the global Needleman-Wunsch alignment. Fill a matrix H where H[i][j] is the score of the best local alignment ending at position i of the first sequence and j of the second:

H[i][j] = max(
    0,                                  // start a new alignment here
    H[i-1][j-1] + s(a_i, b_j),          // match or mismatch
    H[i-1][j]   - gap_penalty,          // gap in sequence b
    H[i][j-1]   - gap_penalty           // gap in sequence a
)

The zero is what makes it local: any prefix that has accumulated a negative score is abandoned rather than carried. The best local alignment is the largest value anywhere in the matrix, and the alignment itself is read out by tracing back from that cell until you reach a zero.

The substitution score s is not a match/mismatch indicator for proteins. It comes from a matrix such as BLOSUM62, whose entries are log-odds ratios of how often a residue pair is observed in aligned homologous blocks against how often it would occur by chance, so a leucine aligned to an isoleucine scores positive and a leucine aligned to a proline scores strongly negative. Gaps are usually affine — one cost to open, a smaller cost per residue to extend — which requires two extra matrices but the same recursion shape.

The problem is the cost: filling the matrix is order m times n cell updates for sequences of length m and n. Against a database, n is the concatenated length of every sequence in it, which is on the order of hundreds of billions of residues for a comprehensive protein database. Running an exact algorithm against that per query is the thing BLAST exists to not do.

A worked scoring matrix

A small case makes the recursion concrete. Align GATTACA and GTTACA with a simple scheme: match plus 2, mismatch minus 1, gap minus 2. Filling the first few rows:

          -   G   T   T   A   C   A
      -   0   0   0   0   0   0   0
      G   0   2   0   0   0   0   0
      A   0   0   1   0   2   0   2
      T   0   0   2   3   1   1   0
      T   0   0   2   4   2   0   0
      A   0   0   0   2   6   4   2
      C   0   0   0   0   4   8   6
      A   0   0   0   0   2   6  10

The maximum is 10 in the bottom-right cell. Tracing back diagonally from it recovers ATTACA aligned to -TTACA with the leading G matched and one deletion: the six matching residues score 12 and the single gap costs 2. The mechanism to notice is that the score builds monotonically along a diagonal of matches and collapses toward zero away from one. That behaviour — high-scoring regions live on diagonals — is exactly the observation the heuristic exploits.

Seed and extend

BLAST, published by Stephen Altschul and colleagues in the Journal of Molecular Biology in 1990, replaces the exhaustive matrix fill with three phases.

1. Seeding. Break the query into every substring of length W, the word size. For each word, find the database positions where it occurs. For nucleotides this is exact word matching. For proteins it is not: BLAST generates, for each query word, the list of all words scoring at least a threshold T against it under the substitution matrix, and looks up those too. So a query word LKV also seeds on IKV and LRV, because those score highly in BLOSUM62. This neighbourhood step is why protein BLAST finds distant homologues that exact word matching would miss, and it is the part usually omitted from a summary of how BLAST works.
2. Extension. A seed hit is extended outward in both directions without gaps, accumulating score, and abandoned when the score drops more than X below the best value seen so far. That X-drop cutoff is what keeps extension from running the full length of the database sequence. The result is a high-scoring segment pair.
3. Gapped alignment. Segment pairs that pass a score threshold get a real gapped dynamic-programming alignment, but only in a band around the seed rather than over the whole matrix. The 1997 gapped BLAST and PSI-BLAST paper by Altschul and colleagues in Nucleic Acids Research added a further requirement: two non- overlapping word hits on the same diagonal within a window, before extension is attempted at all. Requiring two hits raises the evidence needed to spend the extension cost and made the program substantially faster at the same sensitivity.

The trade is stated in one sentence: an alignment whose best region contains no word of length W scoring above T is invisible to BLAST no matter how good the alignment would have been. That is the sensitivity you buy the speed with, and it is not a rounding error — it is the reason more sensitive methods exist for remote homology, and the reason embedding-based sequence search is interesting.

What an E-value is a function of

NCBI defines the Expect value as the number of hits scoring at least this well that you would expect to see by chance in a database of this size. The Karlin-Altschul statistics behind it give

E = K * m * n * exp(-lambda * S)

  m       effective query length
  n       effective database length (total residues)
  S       raw alignment score
  K, lambda   constants determined by the scoring matrix
              and the residue composition

Normalising the raw score by those two constants gives the bit score, S-prime equal to lambda times S minus the natural log of K, all divided by the log of 2, which collapses the expression to

E = m * n * 2^(-S')

Work an example. A 300-residue query, a database of 2 times 10 to the eleventh residues, a hit at 60 bits. Then E equals 300 times 2e11 divided by 2 to the sixtieth. Two to the sixtieth is about 1.15 times 10 to the eighteenth, and 300 times 2e11 is 6 times 10 to the thirteenth, so E is about 5 times 10 to the minus 5.

Two consequences fall straight out of that formula and both are routinely missed. First, E is proportional to database size, so the same alignment against a database ten times larger has an E-value ten times worse — the alignment did not change, the multiple-testing burden did. Second, the bit score is not a function of database size at all, which is why bit score is the right quantity to compare across searches and E-value is the right quantity to threshold within one. A hit at E equal to 1e-5 in 2005 may be at E equal to 1e-3 today purely because the databases grew.

The parameters that change your results

• Word size. The defaults differ sharply by program. NCBI’s BLAST+ documentation gives blastp a word size of 3, the megablast task a word size of 28, and the -task blastn setting a word size of 11 — and the bare blastn command defaults to the megablast task, which is the trap. If you run blastn expecting classic behaviour and get few hits on a cross-species search, you were running megablast at word size 28, which is tuned for near-identical sequences. Pass -task blastn for divergent nucleotide comparisons.
• The substitution matrix. BLOSUM62 is the default and is tuned for moderate divergence. For short queries or close homologues a matrix built from more similar blocks is more appropriate, and for remote homology a more permissive one is. The matrix also changes K and lambda, and therefore every E-value.
• Low-complexity filtering. Regions of biased composition generate high-scoring alignments that mean nothing, so they are masked by default. Turning masking off will produce impressive-looking hits between two poly-glutamine stretches.
• The E-value cutoff is not a p-value. It is an expected count, so values above 1 are meaningful (you expect several such hits by chance) in a way a probability above 1 could not be.