Statistics of Sequence Alignment and BLAST Scores

When you align two sequences, you get a numeric score — a sum of substitution rewards/penalties (from a matrix like BLOSUM62) minus gap penalties, computed by dynamic programming for global alignment (Needleman-Wunsch Algorithm) or local alignment (Smith–Waterman). The question that matters biologically is never "what's the score?" but "is this score higher than we'd expect from unrelated sequences?" That requires knowing the null distribution of alignment scores generated by chance alone. Here is the key statistical fact: if you take two long random sequences and slide a local-alignment algorithm across them, the best local alignment score is the maximum of very many (correlated) segment scores, each of which — for scoring systems with negative expected score per random match — behaves like a rare, exponentially unlikely upward excursion. The theory of extremes says that maxima of many independent (or weakly dependent), identically distributed random variables converge to one of a few universal limiting families as the number of trials grows. For local alignment scores under the conditions Karlin and Altschul specified (negative expected score per aligned pair, positive probability of a positive score), that limiting family is the Gumbel distribution, a particular type of extreme-value distribution with a characteristic long right tail. This is the same mathematical family that describes record-breaking flood heights or maximum daily temperatures — situations where you care about the biggest of many trials, not the typical one. The practical output of this theory has two ingredients: a bit score, which rescales the raw alignment score into matrix-independent, comparable units, and an E-value, which converts the bit score into an expected count of chance hits given the sizes of the sequences being searched. Together they let a biologist ask, precisely, "how many alignments this good would I expect to see by pure chance in a database this size?" — turning an alignment score into an actionable statistical judgment about homology, a concept closely tied to Homology and Sequence Conservation and to the broader logic of Hypothesis Testing in Biology.

Step 1 — Scoring a local alignment. For two sequences, every possible local ungapped alignment (matching contiguous or gapped stretches) is assigned a score using a substitution matrix such as BLOSUM62 for proteins or a match/mismatch scheme for DNA, with gap-open and gap-extension penalties subtracted for insertions/deletions. Smith–Waterman dynamic programming finds the optimal (maximum-scoring) local alignment exactly; BLAST uses a fast heuristic (seed-and-extend via short exact-match "words") that approximates this optimum without exhaustively searching all alignments, trading a little sensitivity for enormous speed.

Step 2 — Why chance scores follow an extreme-value law. Consider comparing two random sequences of length m and n under a scoring scheme where the expected score of a randomly aligned pair of residues is negative (essential — otherwise arbitrarily long alignments would always improve the score, and there would be no meaningful "best local score"). Under this condition, local alignment scores computed at many different starting points behave like a large collection of relatively independent trials, each usually contributing a low or negative score, but occasionally producing a high positive excursion by chance. The single highest score among a huge number of such trials is a sample maximum, and extreme-value theory tells us the distribution of this maximum converges (as m, n → ∞) to:

This is the Gumbel-type extreme-value distribution. Two parameters, λ (a scale parameter that depends on the substitution matrix and background residue frequencies) and K (a constant related to the effective search space), fully characterize the null model for a given scoring scheme. This result — the Karlin–Altschul equation — was derived rigorously by Samuel Karlin and Stephen Altschul (1990) for ungapped local alignments, and it is the mathematical engine behind every E-value BLAST reports.

Step 3 — Normalizing to bit scores. Raw scores S depend on the arbitrary scale of the substitution matrix (BLOSUM62 scores are integers with a particular scaling). To compare results across different matrices and searches, the raw score is converted to a matrix-independent bit score:

Bit scores are in units of information (bits), which is why two completely different scoring schemes that both correctly reflect the same underlying biology give comparable bit scores — the λ and K normalization absorbs the arbitrary scale.

Step 4 — Converting to an E-value. Once we have the bit score, the expected number of distinct local alignments with score ≥ S′ occurring purely by chance in a search of a query of effective length m′ against a database of effective size n′ is:

This is the E-value. It is not a probability (it can exceed 1); it is an expected count. An E-value of 3 means that, by chance alone, you'd expect about three unrelated sequences in a database this size to score this well or better. An E-value of 10⁻⁴⁵ means you'd need to search an astronomically larger database before chance alone would be expected to produce even one hit that good — strong evidence of true homology.

Step 5 — The P-value connection. For a single alignment (not the best-of-many-comparisons search context), the probability of observing at least one chance alignment with score ≥ S is related to the E-value by the Poisson approximation:

For small E (E < 0.01), P ≈ E, so the two are nearly interchangeable in the regime that matters for declaring significance; for large E they diverge sharply, which is one reason BLAST reports E-values rather than P-values as its primary statistic — E-values remain interpretable (as expected counts) even when not small.

Step 6 — Gapped alignments: theory meets heuristics. The clean Karlin–Altschul theory above is proven rigorously only for ungapped local alignments. Real BLAST alignments include gaps (insertions/deletions), for which no closed-form extreme-value theory exists. Altschul and colleagues showed empirically, through massive simulation with random sequences under realistic gap penalty schemes, that gapped alignment scores also closely follow the same Gumbel-type extreme-value form, just with λ and K parameters estimated empirically (via simulation) for each combination of substitution matrix and gap penalties, rather than derived analytically. This empirical extension — validated by curve-fitting millions of simulated random-sequence alignments — is what actually powers modern BLASTP/BLASTN/BLASTX E-values.

The statistical theory of local alignment scores was not guessed — it was proven mathematically for a restricted case and then validated empirically for the general case actually used in practice.

Karlin and Altschul (1990), Proceedings of the National Academy of Sciences. Samuel Karlin (a probabilist) and Stephen Altschul (then at NIH, later a lead author of BLAST) proved that for ungapped local alignments between two random sequences — generated from independent, identically distributed letters with a scoring scheme whose expected score per aligned pair is negative but which has some positive-scoring outcomes — the distribution of the maximal segment score converges, as sequence lengths grow, to the Gumbel extreme-value form:

They derived closed-form (or numerically computable) expressions for λ and K directly from the scoring matrix and the background amino-acid or nucleotide frequencies used. This was a rigorous asymptotic proof, using tools from random walk theory and large-deviation statistics, not a simulation fit — it showed why the extreme-value law, and not the normal law suggested by naive "sum of many random variables" intuition, is the correct null distribution.

Altschul, Gish, Miller, Myers, and Lipman (1990), and the original BLAST paper. The team that built BLAST used the Karlin–Altschul ungapped theory directly, because the original BLAST algorithm (1990) searched for high-scoring ungapped segment pairs (HSPs). This let them attach rigorous E-values to every reported hit from day one, a major advance over earlier alignment tools like FASTA, which lacked a comparably principled statistical foundation.

Extending to gapped alignments — simulation as evidence. When BLAST2 and later versions incorporated gapped alignments (to better detect homologs with insertions/deletions), no analytic proof of the extreme-value limit existed for gapped scores. Altschul and colleagues (notably with Wolfgang Bucher, and later systematized by Sean Eddy and others) addressed this empirically: they generated enormous numbers of random sequence pairs (with realistic amino-acid composition, drawn to mimic real proteins but scrambled to destroy any true homology), computed optimal gapped local alignment scores for each pair using Smith–Waterman-type dynamic programming with specific gap-open/extend penalties, and fit the resulting empirical score distributions. The consistent finding, replicated across many substitution matrices (BLOSUM45, 62, 80), many gap penalty combinations, and both protein and nucleotide alphabets: gapped local alignment scores from random sequences still closely follow the Gumbel extreme-value distribution, just with λ and K values that must be estimated by curve-fitting to simulated data rather than derived in closed form. This finding — essentially a large-scale empirical hypothesis test of "does the EVD model fit gapped scores as well as it fits ungapped scores?" — was the evidence that let NCBI publish precomputed λ/K lookup tables for standard matrix/gap-penalty combinations, which is what your BLAST search silently uses today.

Corroborating evidence from genome-scale searches. A further, ongoing line of evidence comes from real-world validation: when BLAST is run against genuinely randomized ("shuffled") sequence databases that preserve amino-acid composition but destroy all real biological signal, the empirical distribution of best-hit scores matches the Gumbel prediction closely, and the reported E-values correctly predict the number of false hits observed. This calibration check — repeated whenever new scoring schemes or statistical models are introduced (e.g., for PSI-BLAST's position-specific matrices) — is the ongoing empirical proof that the model keeps working as databases and algorithms evolve.

Example 1 — Converting a raw score to bit score and E-value. Suppose a BLASTP search of a 250-residue query against a protein database of 2 × 10⁹ total residues (a moderate-sized non-redundant database) returns a local alignment with:

• Raw score S = 180 (BLOSUM62 units)
• λ = 0.267, K = 0.041 (typical BLOSUM62/standard-gap-penalty values)
• Effective query length m′ ≈ 220 (query length minus edge-effect correction)
• Effective database size n′ ≈ 2 × 10⁹

Step 1 — bit score:

Step 2 — E-value:

$2^{74} \approx 1.89 \times 10^{22}$, so:

Interpretation: an E-value of 2.3 × 10⁻¹¹ means that in a database this size, you'd expect roughly 1 chance alignment this good for every 4 × 10¹⁰ searches like this one — vanishingly unlikely to arise without true homology. This is a textbook "highly significant hit," the kind reported for a clear ortholog.

Example 2 — Same bit score, different database: why E-values are not fixed properties of an alignment. Now imagine the identical alignment (same bit score, 74.0 bits) is instead reported from a search against a small, specialized database — say, all 1,200 annotated proteins of a single virus family, totaling n′ ≈ 4 × 10⁵ residues.

The E-value becomes even smaller (more significant) — because with fewer sequences searched, a hit this strong is even less likely to be a coincidence. Conversely, if the same alignment came from searching the entire 3 × 10¹¹-residue non-redundant "nr" database (much larger than Example 1's hypothetical), the E-value would inflate roughly proportionally to the database size increase — a bit score of 74 might correspond to E ≈ 10⁻⁹ rather than 10⁻¹¹.

The lesson: the bit score is a property of the alignment itself (matrix-independent, comparable across searches). The E-value is a property of the alignment in the context of a specific search — it depends on database size and query length. This is why researchers comparing hits across different BLAST runs, or across different versions of a growing database, should compare bit scores, not raw E-values, and why a "good E-value" from a niche database search does not carry over unchanged to a genome-wide search.

Samuel Karlin (1924–2007) was a mathematician at Stanford, renowned for foundational work in probability theory, mathematical population genetics, and — later in his career — molecular sequence statistics. Trained originally in pure and applied mathematics, Karlin brought rigorous large-deviation and random-walk theory to bear on a problem biologists had been solving heuristically: what does a "good" alignment score actually mean statistically? Working with Stephen Altschul, then a young statistician-biologist at the NIH who would go on to co-author BLAST itself, Karlin supplied the theorems that turned alignment scoring from an ad hoc heuristic into a calibrated statistical test.

Genome annotation pipelines run BLAST (or its faster descendants, DIAMOND, MMseqs2) against every predicted gene in a newly sequenced genome, using E-value cutoffs (commonly 10⁻⁵ or 10⁻¹⁰) to decide which genes get functional labels transferred from known homologs — a decision repeated for every genome deposited in GenBank, RefSeq, and UniProt.

Clinical microbiology and outbreak investigation use BLAST-style alignment statistics to identify pathogens from sequencing reads, matching a patient's isolate against reference databases and using E-values to distinguish a genuine species/strain match from a spurious low-complexity match — critical during rapid outbreak response (e.g., identifying a novel coronavirus lineage from metagenomic sequencing).

Drug target and resistance-gene discovery relies on E-value-filtered BLAST hits to flag candidate antibiotic-resistance genes or virulence factors in bacterial genomes, feeding directly into surveillance systems like CARD (Comprehensive Antibiotic Resistance Database).

Protein structure and function prediction tools (PSI-BLAST, HMMER, and profile-based successors) build on the same extreme-value statistical framework, iteratively refining position-specific scoring matrices while still reporting Karlin–Altschul-style E-values at each round — a technique essential for finding remote homologs whose sequence identity has decayed below the reach of a single BLAST pass.

Conservation genomics and biodiversity surveys use BLAST against curated reference databases (e.g., barcoding loci like COI or 16S rRNA) to assign species identity from environmental DNA samples, with E-value thresholds tuned to the reference database's size and quality — directly informing biodiversity monitoring and invasive-species detection.

Foreshadowing: the same logic of "is this statistic more extreme than chance would produce across many simultaneous comparisons" reappears, in a different guise, in Genome-Wide Association Statistics, where millions of simultaneous SNP-association tests demand corrections (genomic control, Bonferroni, FDR) analogous in spirit to the multiple-comparison reasoning embedded in BLAST's E-value — both fields confront the statistics of picking out real signal from the extreme tail of many chance comparisons.

Ask the tutor: How do E-value thresholds need to be adjusted when comparing short reads (e.g., 100 bp sequencing reads) against a reference genome, versus comparing full-length protein queries against a protein database?