BLAST

Introduction

You have just sequenced a stretch of DNA and want to know what it is. Has anyone seen something like it before? In what organism? Does it resemble a known gene? The answer lies in public databases containing many billions of letters, and you want it in seconds.

The textbook tool for comparing two sequences is the Smith–Waterman algorithm (1981), a dynamic program that finds the best local alignment exactly. It fills a table with one cell for every pair of positions. For a query of $m$ letters and a database of $n$ letters that is $m \cdot n$ cells, every time. Doing that for millions of queries against a growing database is out of reach.

BLAST (Basic Local Alignment Search Tool), published in 1990 by Stephen Altschul, Warren Gish, Webb Miller, Eugene Myers and David Lipman, made a bet: any interesting similarity almost certainly contains a short stretch that matches exactly. So look for those short stretches first, and spend effort only around them.

This is the seed-and-extend strategy. It is a close relative of the idea in string matching: use an index to jump straight to candidate positions instead of testing all of them. The price is that BLAST is a heuristic — it is not guaranteed to find the optimal alignment. The reward is a search that is orders of magnitude cheaper, with statistics that tell you how much to trust each hit.

Here is the whole algorithm in pseudocode:

BLAST(query, database, w, X):
  index every w-letter word of the query
  for each sequence S in the database:
    for each position j in S:
      if the w-letter word at j occurs in the query:   # a seed hit
        extend the match left and right without gaps,
          stopping when the score falls X below its best value
        keep the segment if its score is high enough    # an HSP
  report the HSPs ranked by E-value

Search a Database

Below, an 80-letter query is searched against eight random sequences. Three of them (marked ★) secretly contain a copy of the query with about 8%, 20% and 30% of the letters mutated, plus the occasional insertion or deletion.

<p class="hint">A random 80-letter DNA query is searched against 8 random sequences of 300 letters. Three of them secretly contain a mutated copy of the query (★). Change the word size and watch speed and sensitivity trade off.</p>
<div class="ctl">
  <label for="w">Word size w: <b id="wv">8</b></label>
  <input id="w" type="range" min="4" max="12" value="8" aria-label="Word size w">
</div>
<div class="btns">
  <button id="p11" type="button">Fast (w=11)</button>
  <button id="p8" type="button">Balanced (w=8)</button>
  <button id="p5" type="button">Sensitive (w=5)</button>
  <button id="newdata" type="button" class="ghost">New random data</button>
</div>
<div class="cap" id="cap"></div>
<canvas id="plot" width="600" height="160" aria-label="Dot plot of seed hits and extended alignments between the query and the selected sequence"></canvas>
<div class="legend"><span class="k s"></span>seed hit <span class="k h"></span>extended HSP</div>
<div id="sum" class="sum"></div>
<table id="tbl">
  <thead><tr><th>Seq</th><th>Seeds</th><th>HSP score</th><th>E-value</th><th>SW score</th></tr></thead>
  <tbody></tbody>
</table>
<p class="note">★ marks a sequence that really contains a mutated copy of the query. E = K·m·n·e^(−λS) with K = 0.46 and λ computed for these scores. The SW column is the exact optimum, computed for comparison.</p>

:root { --fg:#222; --mut:#555; --bg:#fff; --card:#f1f4f8; --line:#d3dae3; --pri:#1d3557; --seed:#9aa7b6; --hsp:#e63946; --ok:#0a7d33; }
@media (prefers-color-scheme: dark) {
  :root { --fg:#e6e9ee; --mut:#a9b3c0; --bg:#14181d; --card:#1e252d; --line:#38424e; --pri:#7fa8d9; --seed:#6d7b8c; --hsp:#ff6b78; --ok:#4cc27a; }
}
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: var(--fg); background: var(--bg); margin: 0; }
.hint { font-size: .88rem; color: var(--mut); margin: 0 0 .6rem; line-height: 1.45; }
.ctl { display: flex; flex-direction: column; gap: .2rem; margin-bottom: .5rem; font-size: .9rem; }
.ctl input { width: 100%; }
.btns { display: flex; gap: .4rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .7rem; border: 1px solid var(--pri); background: var(--pri); color: var(--bg); border-radius: 8px; cursor: pointer; }
button.ghost { background: transparent; color: var(--pri); }
.cap { font-size: .8rem; color: var(--mut); margin: .2rem 0; }
canvas { width: 100%; height: auto; max-width: 600px; display: block; background: var(--card); border: 1px solid var(--line); border-radius: 6px; }
.legend { font-size: .78rem; color: var(--mut); margin: .3rem 0; }
.k { display: inline-block; width: 10px; height: 10px; margin: 0 4px 0 8px; border-radius: 2px; }
.k.s { background: var(--seed); }
.k.h { background: var(--hsp); }
.sum { font-size: .85rem; line-height: 1.5; background: var(--card); border-radius: 8px; padding: .5rem .7rem; margin: .4rem 0; }
.sum b { color: var(--pri); }
table { border-collapse: collapse; width: 100%; font-size: .78rem; }
th, td { padding: .25rem .3rem; text-align: right; border-bottom: 1px solid var(--line); }
th:first-child, td:first-child { text-align: left; }
tbody tr { cursor: pointer; }
tbody tr.sel { background: var(--card); font-weight: 600; }
td.ok { color: var(--ok); }
.note { font-size: .78rem; color: var(--mut); line-height: 1.4; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Scoring and sizes: +1 match, -2 mismatch, X-drop 10
const MATCH = 1, MISMATCH = -2, GAP = -4, XDROP = 10;
const QLEN = 80, SLEN = 300, NSUB = 8, KA_K = 0.46, E_CUT = 0.001;
const IDENT = [0.92, 0.80, 0.70];
// Small seeded random generator (reproducible data)
function rng(seed) {
  let a = seed >>> 0;
  return () => {
    a = (a + 0x6D2B79F5) >>> 0;
    let t = a;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}
// Karlin-Altschul lambda: the positive root of sum p_i p_j exp(lambda s_ij) = 1
function karlinLambda() {
  let lo = 0.1, hi = 5;
  for (let it = 0; it < 80; it++) {
    const mid = (lo + hi) / 2;
    const f = 0.25 * Math.exp(MATCH * mid) + 0.75 * Math.exp(MISMATCH * mid) - 1;
    if (f > 0) hi = mid; else lo = mid;
  }
  return (lo + hi) / 2;
}
const LAMBDA = karlinLambda();
// Random database with three planted, mutated copies of the query
function makeData(seed) {
  const r = rng(seed), A = 'ACGT';
  const rnd = n => { let s = ''; for (let i = 0; i < n; i++) s += A[Math.floor(r() * 4)]; return s; };
  const q = rnd(QLEN);
  const subs = [];
  for (let i = 0; i < NSUB; i++) subs.push(rnd(SLEN));
  const order = [...Array(NSUB).keys()];
  for (let i = NSUB - 1; i > 0; i--) { const j = Math.floor(r() * (i + 1)); [order[i], order[j]] = [order[j], order[i]]; }
  const planted = [];
  IDENT.forEach((id, n) => {
    let h = '';
    for (let i = 0; i < q.length; i++) {
      const x = r();
      if (x < 0.008) continue;
      if (x < 0.016) h += A[Math.floor(r() * 4)];
      if (r() < id) h += q[i]; else h += A[(A.indexOf(q[i]) + 1 + Math.floor(r() * 3)) % 4];
    }
    const si = order[n], s = subs[si], pos = 20 + Math.floor(r() * (SLEN - h.length - 40));
    subs[si] = s.slice(0, pos) + h + s.slice(pos + h.length);
    planted.push({ si, ident: id });
  });
  planted.sort((a, b) => a.si - b.si);
  return { q, subs, planted };
}
// BLAST-style search: seed words, then ungapped X-drop extension
function blast(q, subs, w) {
  const idx = new Map();
  for (let i = 0; i + w <= q.length; i++) {
    const k = q.substr(i, w);
    if (!idx.has(k)) idx.set(k, []);
    idx.get(k).push(i);
  }
  const N = subs.reduce((a, s) => a + s.length, 0);
  let seeds = 0, ext = 0, cells = 0;
  const per = subs.map(s => {
    const lastEnd = new Map(), hits = [], hsps = [];
    for (let j = 0; j + w <= s.length; j++) {
      const list = idx.get(s.substr(j, w));
      if (!list) continue;
      for (const i of list) {
        seeds++; hits.push([i, j]);
        const d = j - i;
        if ((lastEnd.get(d) ?? -1) > j) continue;
        ext++; cells += w;
        let run = 0, gr = 0, lr = 0;
        for (let a = i + w, b = j + w, n = 1; a < q.length && b < s.length; a++, b++, n++) {
          run += q[a] === s[b] ? MATCH : MISMATCH; cells++;
          if (run > gr) { gr = run; lr = n; }
          if (gr - run > XDROP) break;
        }
        run = 0; let gl = 0, ll = 0;
        for (let a = i - 1, b = j - 1, n = 1; a >= 0 && b >= 0; a--, b--, n++) {
          run += q[a] === s[b] ? MATCH : MISMATCH; cells++;
          if (run > gl) { gl = run; ll = n; }
          if (gl - run > XDROP) break;
        }
        lastEnd.set(d, j + w + lr);
        hsps.push({ qs: i - ll, ss: j - ll, len: w + lr + ll, score: w * MATCH + gr + gl });
      }
    }
    hsps.sort((x, y) => y.score - x.score);
    const best = hsps[0] || null;
    const E = best ? KA_K * q.length * N * Math.exp(-LAMBDA * best.score) : Infinity;
    return { hits, hsps, best, E, seeds: hits.length };
  });
  return { per, seeds, ext, cells, N };
}
// Smith-Waterman baseline: every cell of the table is filled
function swScore(q, s) {
  let prev = new Int32Array(s.length + 1), cur = new Int32Array(s.length + 1), best = 0;
  for (let i = 1; i <= q.length; i++) {
    for (let j = 1; j <= s.length; j++) {
      const d = prev[j - 1] + (q[i - 1] === s[j - 1] ? MATCH : MISMATCH);
      const v = Math.max(0, d, prev[j] + GAP, cur[j - 1] + GAP);
      cur[j] = v; if (v > best) best = v;
    }
    [prev, cur] = [cur, prev];
  }
  return best;
}
// Interface (DOM)
const $ = id => document.getElementById(id);
let seed = 7, w = 8, sel = 0, data, res, sw;
const fmtE = e => !isFinite(e) ? '—' : e < 0.01 ? e.toExponential(1) : e.toFixed(2);
const css = n => getComputedStyle(document.documentElement).getPropertyValue(n);
function draw() {
  const c = $('plot'), g = c.getContext('2d'), sx = c.width / SLEN, sy = c.height / QLEN;
  g.clearRect(0, 0, c.width, c.height);
  g.fillStyle = css('--seed');
  res.per[sel].hits.forEach(([i, j]) => g.fillRect(j * sx, i * sy, Math.max(2, sx), Math.max(2, sy)));
  g.strokeStyle = css('--hsp'); g.lineWidth = 3;
  res.per[sel].hsps.filter(h => h.score >= 12).forEach(h => {
    g.beginPath(); g.moveTo(h.ss * sx, h.qs * sy); g.lineTo((h.ss + h.len) * sx, (h.qs + h.len) * sy); g.stroke();
  });
  $('cap').textContent = i18nFmt("Dot plot: sequence %n% (x) against the query (y), word size %w%. Click a row to change sequence.", { n: sel + 1, w });
}
function run() {
  $('wv').textContent = w;
  res = blast(data.q, data.subs, w);
  sw = data.subs.map(s => swScore(data.q, s));
  const full = QLEN * res.N;
  const isP = i => data.planted.some(p => p.si === i);
  const found = data.planted.filter(p => res.per[p.si].E <= E_CUT).length;
  const fa = res.per.filter((p, i) => !isP(i) && p.E <= E_CUT).length;
  $('sum').innerHTML =
    i18nFmt("Seed hits: <b>%s%</b> · extensions started: <b>%e%</b>", { s: res.seeds, e: res.ext }) + '<br>' +
    i18nFmt("Cells computed — BLAST: <b>%b%</b> · full Smith–Waterman: <b>%f%</b> (%x%× fewer)", { b: res.cells.toLocaleString(), f: full.toLocaleString(), x: (full / Math.max(1, res.cells)).toFixed(1) }) + '<br>' +
    i18nFmt("Planted homologs found (E ≤ 0.001): <b>%f%</b> of %t% · false alarms: <b>%fa%</b>", { f: found, t: data.planted.length, fa });
  const tb = document.querySelector('#tbl tbody'); tb.innerHTML = '';
  res.per.forEach((p, i) => {
    const tr = document.createElement('tr');
    if (i === sel) tr.className = 'sel';
    const hit = p.E <= E_CUT;
    tr.innerHTML = '<td>' + (i + 1) + (isP(i) ? ' ★' : '') + '</td><td>' + p.seeds + '</td><td>' + (p.best ? p.best.score : '—') +
      '</td><td class="' + (hit ? 'ok' : '') + '">' + fmtE(p.E) + '</td><td>' + sw[i] + '</td>';
    tr.addEventListener('click', () => { sel = i; run(); });
    tb.appendChild(tr);
  });
  draw();
}
function fresh() { data = makeData(seed); sel = data.planted[0].si; run(); }
$('w').addEventListener('input', e => { w = +e.target.value; run(); });
const preset = v => () => { w = v; $('w').value = v; run(); };
$('p11').addEventListener('click', preset(11));
$('p8').addEventListener('click', preset(8));
$('p5').addEventListener('click', preset(5));
$('newdata').addEventListener('click', () => { seed = seed * 7 + 3; fresh(); });
fresh();

<p class="hint">A random 80-letter DNA query is searched against 8 random sequences of 300 letters. Three of them secretly contain a mutated copy of the query (★). Change the word size and watch speed and sensitivity trade off.</p>
<div class="ctl">
  <label for="w">Word size w: <b id="wv">8</b></label>
  <input id="w" type="range" min="4" max="12" value="8" aria-label="Word size w">
</div>
<div class="btns">
  <button id="p11" type="button">Fast (w=11)</button>
  <button id="p8" type="button">Balanced (w=8)</button>
  <button id="p5" type="button">Sensitive (w=5)</button>
  <button id="newdata" type="button" class="ghost">New random data</button>
</div>
<div class="cap" id="cap"></div>
<canvas id="plot" width="600" height="160" aria-label="Dot plot of seed hits and extended alignments between the query and the selected sequence"></canvas>
<div class="legend"><span class="k s"></span>seed hit <span class="k h"></span>extended HSP</div>
<div id="sum" class="sum"></div>
<table id="tbl">
  <thead><tr><th>Seq</th><th>Seeds</th><th>HSP score</th><th>E-value</th><th>SW score</th></tr></thead>
  <tbody></tbody>
</table>
<p class="note">★ marks a sequence that really contains a mutated copy of the query. E = K·m·n·e^(−λS) with K = 0.46 and λ computed for these scores. The SW column is the exact optimum, computed for comparison.</p>

:root { --fg:#222; --mut:#555; --bg:#fff; --card:#f1f4f8; --line:#d3dae3; --pri:#1d3557; --seed:#9aa7b6; --hsp:#e63946; --ok:#0a7d33; }
@media (prefers-color-scheme: dark) {
  :root { --fg:#e6e9ee; --mut:#a9b3c0; --bg:#14181d; --card:#1e252d; --line:#38424e; --pri:#7fa8d9; --seed:#6d7b8c; --hsp:#ff6b78; --ok:#4cc27a; }
}
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: var(--fg); background: var(--bg); margin: 0; }
.hint { font-size: .88rem; color: var(--mut); margin: 0 0 .6rem; line-height: 1.45; }
.ctl { display: flex; flex-direction: column; gap: .2rem; margin-bottom: .5rem; font-size: .9rem; }
.ctl input { width: 100%; }
.btns { display: flex; gap: .4rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .7rem; border: 1px solid var(--pri); background: var(--pri); color: var(--bg); border-radius: 8px; cursor: pointer; }
button.ghost { background: transparent; color: var(--pri); }
.cap { font-size: .8rem; color: var(--mut); margin: .2rem 0; }
canvas { width: 100%; height: auto; max-width: 600px; display: block; background: var(--card); border: 1px solid var(--line); border-radius: 6px; }
.legend { font-size: .78rem; color: var(--mut); margin: .3rem 0; }
.k { display: inline-block; width: 10px; height: 10px; margin: 0 4px 0 8px; border-radius: 2px; }
.k.s { background: var(--seed); }
.k.h { background: var(--hsp); }
.sum { font-size: .85rem; line-height: 1.5; background: var(--card); border-radius: 8px; padding: .5rem .7rem; margin: .4rem 0; }
.sum b { color: var(--pri); }
table { border-collapse: collapse; width: 100%; font-size: .78rem; }
th, td { padding: .25rem .3rem; text-align: right; border-bottom: 1px solid var(--line); }
th:first-child, td:first-child { text-align: left; }
tbody tr { cursor: pointer; }
tbody tr.sel { background: var(--card); font-weight: 600; }
td.ok { color: var(--ok); }
.note { font-size: .78rem; color: var(--mut); line-height: 1.4; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Scoring and sizes: +1 match, -2 mismatch, X-drop 10
const MATCH = 1, MISMATCH = -2, GAP = -4, XDROP = 10;
const QLEN = 80, SLEN = 300, NSUB = 8, KA_K = 0.46, E_CUT = 0.001;
const IDENT = [0.92, 0.80, 0.70];
// Small seeded random generator (reproducible data)
function rng(seed) {
  let a = seed >>> 0;
  return () => {
    a = (a + 0x6D2B79F5) >>> 0;
    let t = a;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}
// Karlin-Altschul lambda: the positive root of sum p_i p_j exp(lambda s_ij) = 1
function karlinLambda() {
  let lo = 0.1, hi = 5;
  for (let it = 0; it < 80; it++) {
    const mid = (lo + hi) / 2;
    const f = 0.25 * Math.exp(MATCH * mid) + 0.75 * Math.exp(MISMATCH * mid) - 1;
    if (f > 0) hi = mid; else lo = mid;
  }
  return (lo + hi) / 2;
}
const LAMBDA = karlinLambda();
// Random database with three planted, mutated copies of the query
function makeData(seed) {
  const r = rng(seed), A = 'ACGT';
  const rnd = n => { let s = ''; for (let i = 0; i < n; i++) s += A[Math.floor(r() * 4)]; return s; };
  const q = rnd(QLEN);
  const subs = [];
  for (let i = 0; i < NSUB; i++) subs.push(rnd(SLEN));
  const order = [...Array(NSUB).keys()];
  for (let i = NSUB - 1; i > 0; i--) { const j = Math.floor(r() * (i + 1)); [order[i], order[j]] = [order[j], order[i]]; }
  const planted = [];
  IDENT.forEach((id, n) => {
    let h = '';
    for (let i = 0; i < q.length; i++) {
      const x = r();
      if (x < 0.008) continue;
      if (x < 0.016) h += A[Math.floor(r() * 4)];
      if (r() < id) h += q[i]; else h += A[(A.indexOf(q[i]) + 1 + Math.floor(r() * 3)) % 4];
    }
    const si = order[n], s = subs[si], pos = 20 + Math.floor(r() * (SLEN - h.length - 40));
    subs[si] = s.slice(0, pos) + h + s.slice(pos + h.length);
    planted.push({ si, ident: id });
  });
  planted.sort((a, b) => a.si - b.si);
  return { q, subs, planted };
}
// BLAST-style search: seed words, then ungapped X-drop extension
function blast(q, subs, w) {
  const idx = new Map();
  for (let i = 0; i + w <= q.length; i++) {
    const k = q.substr(i, w);
    if (!idx.has(k)) idx.set(k, []);
    idx.get(k).push(i);
  }
  const N = subs.reduce((a, s) => a + s.length, 0);
  let seeds = 0, ext = 0, cells = 0;
  const per = subs.map(s => {
    const lastEnd = new Map(), hits = [], hsps = [];
    for (let j = 0; j + w <= s.length; j++) {
      const list = idx.get(s.substr(j, w));
      if (!list) continue;
      for (const i of list) {
        seeds++; hits.push([i, j]);
        const d = j - i;
        if ((lastEnd.get(d) ?? -1) > j) continue;
        ext++; cells += w;
        let run = 0, gr = 0, lr = 0;
        for (let a = i + w, b = j + w, n = 1; a < q.length && b < s.length; a++, b++, n++) {
          run += q[a] === s[b] ? MATCH : MISMATCH; cells++;
          if (run > gr) { gr = run; lr = n; }
          if (gr - run > XDROP) break;
        }
        run = 0; let gl = 0, ll = 0;
        for (let a = i - 1, b = j - 1, n = 1; a >= 0 && b >= 0; a--, b--, n++) {
          run += q[a] === s[b] ? MATCH : MISMATCH; cells++;
          if (run > gl) { gl = run; ll = n; }
          if (gl - run > XDROP) break;
        }
        lastEnd.set(d, j + w + lr);
        hsps.push({ qs: i - ll, ss: j - ll, len: w + lr + ll, score: w * MATCH + gr + gl });
      }
    }
    hsps.sort((x, y) => y.score - x.score);
    const best = hsps[0] || null;
    const E = best ? KA_K * q.length * N * Math.exp(-LAMBDA * best.score) : Infinity;
    return { hits, hsps, best, E, seeds: hits.length };
  });
  return { per, seeds, ext, cells, N };
}
// Smith-Waterman baseline: every cell of the table is filled
function swScore(q, s) {
  let prev = new Int32Array(s.length + 1), cur = new Int32Array(s.length + 1), best = 0;
  for (let i = 1; i <= q.length; i++) {
    for (let j = 1; j <= s.length; j++) {
      const d = prev[j - 1] + (q[i - 1] === s[j - 1] ? MATCH : MISMATCH);
      const v = Math.max(0, d, prev[j] + GAP, cur[j - 1] + GAP);
      cur[j] = v; if (v > best) best = v;
    }
    [prev, cur] = [cur, prev];
  }
  return best;
}
// Interface (DOM)
const $ = id => document.getElementById(id);
let seed = 7, w = 8, sel = 0, data, res, sw;
const fmtE = e => !isFinite(e) ? '—' : e < 0.01 ? e.toExponential(1) : e.toFixed(2);
const css = n => getComputedStyle(document.documentElement).getPropertyValue(n);
function draw() {
  const c = $('plot'), g = c.getContext('2d'), sx = c.width / SLEN, sy = c.height / QLEN;
  g.clearRect(0, 0, c.width, c.height);
  g.fillStyle = css('--seed');
  res.per[sel].hits.forEach(([i, j]) => g.fillRect(j * sx, i * sy, Math.max(2, sx), Math.max(2, sy)));
  g.strokeStyle = css('--hsp'); g.lineWidth = 3;
  res.per[sel].hsps.filter(h => h.score >= 12).forEach(h => {
    g.beginPath(); g.moveTo(h.ss * sx, h.qs * sy); g.lineTo((h.ss + h.len) * sx, (h.qs + h.len) * sy); g.stroke();
  });
  $('cap').textContent = i18nFmt("Dot plot: sequence %n% (x) against the query (y), word size %w%. Click a row to change sequence.", { n: sel + 1, w });
}
function run() {
  $('wv').textContent = w;
  res = blast(data.q, data.subs, w);
  sw = data.subs.map(s => swScore(data.q, s));
  const full = QLEN * res.N;
  const isP = i => data.planted.some(p => p.si === i);
  const found = data.planted.filter(p => res.per[p.si].E <= E_CUT).length;
  const fa = res.per.filter((p, i) => !isP(i) && p.E <= E_CUT).length;
  $('sum').innerHTML =
    i18nFmt("Seed hits: <b>%s%</b> · extensions started: <b>%e%</b>", { s: res.seeds, e: res.ext }) + '<br>' +
    i18nFmt("Cells computed — BLAST: <b>%b%</b> · full Smith–Waterman: <b>%f%</b> (%x%× fewer)", { b: res.cells.toLocaleString(), f: full.toLocaleString(), x: (full / Math.max(1, res.cells)).toFixed(1) }) + '<br>' +
    i18nFmt("Planted homologs found (E ≤ 0.001): <b>%f%</b> of %t% · false alarms: <b>%fa%</b>", { f: found, t: data.planted.length, fa });
  const tb = document.querySelector('#tbl tbody'); tb.innerHTML = '';
  res.per.forEach((p, i) => {
    const tr = document.createElement('tr');
    if (i === sel) tr.className = 'sel';
    const hit = p.E <= E_CUT;
    tr.innerHTML = '<td>' + (i + 1) + (isP(i) ? ' ★' : '') + '</td><td>' + p.seeds + '</td><td>' + (p.best ? p.best.score : '—') +
      '</td><td class="' + (hit ? 'ok' : '') + '">' + fmtE(p.E) + '</td><td>' + sw[i] + '</td>';
    tr.addEventListener('click', () => { sel = i; run(); });
    tb.appendChild(tr);
  });
  draw();
}
function fresh() { data = makeData(seed); sel = data.planted[0].si; run(); }
$('w').addEventListener('input', e => { w = +e.target.value; run(); });
const preset = v => () => { w = v; $('w').value = v; run(); };
$('p11').addEventListener('click', preset(11));
$('p8').addEventListener('click', preset(8));
$('p5').addEventListener('click', preset(5));
$('newdata').addEventListener('click', () => { seed = seed * 7 + 3; fresh(); });
fresh();

Result

Console

Start with Balanced (w=8). The dot plot shows every seed hit as a grey dot. Most are chance matches scattered at random, but a real homolog lights up a diagonal of seeds, and the red line is the extended segment (the HSP, high-scoring segment pair). Click a row to inspect another sequence.

Now choose Fast (w=11): the number of computed cells drops, but the most diverged homolog starts to slip away because it no longer contains an 11-letter exact match. Choose Sensitive (w=5): the weakest copy is found, but there are many more chance seeds, so the work grows. Notice what the E-value column does: the planted copies get tiny E-values, while the random sequences stay far above the threshold, even when they collect many chance seeds.

The Real Complexity

The exact baseline. Smith–Waterman takes $\Theta(mn)$ time in the best, average and worst case, because it fills the whole table. It needs only $O(\min(m,n))$ memory if you want the score alone, or $O(m+n)$ with Hirschberg's trick if you also want the alignment. Comparing against a database of total length $N$ costs $\Theta(mN)$ .

BLAST on random DNA. Take the four letters as uniform and independent. Indexing the query costs $O(m)$ and scanning the database costs $O(N)$ . A given pair of positions forms a seed hit with probability $4^{-w}$ , so the expected number of seed hits is

\mathbb{E}[\text{hits}] = (m-w+1)(N-w+1)\,4^{-w} \approx \frac{mN}{4^{w}}.

Each hit triggers an extension whose expected length is constant: with a match score of $+1$ and a mismatch of $-2$ a random diagonal drifts downward, and the X-drop rule stops it quickly. The expected total work is therefore about $O\!\left(m + N + \dfrac{mN}{4^{w}}\right)$ . Every extra letter of $w$ divides the extension work by 4, and the demo's cell counter shows exactly this.

The worst case is still $O(mN)$ : a repetitive query (say, a long run of one letter) makes almost every position a seed. That is why real programs mask low-complexity regions before searching.

Sensitivity is a probability, not a guarantee. If a homolog has identity $p$ per position, a given $w$ -letter window is an exact match with probability $p^{w}$ . A region of length $L$ has $L-w+1$ windows, so the chance of at least one seed is roughly $1-(1-p^{w})^{L-w+1}$ (an approximation, since neighbouring windows are correlated). With $p=0.7$ and $w=11$ each window succeeds with probability only $0.7^{11}\approx 2\%$ , and a weak homolog can be missed entirely. With $w=5$ it succeeds with $0.7^{5}\approx 17\%$ , at the price of $4^{6}=4096$ times more chance hits.

E-values. Karlin and Altschul (1990) showed that for local ungapped alignments with negative expected score, the best score between random sequences of lengths $m$ and $n$ follows an extreme-value law. The expected number of segments scoring at least $S$ by chance is

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

where $\lambda$ is the unique positive solution of $\sum_{i,j} p_i p_j\, e^{\lambda s_{ij}} = 1$ and $K$ is a constant depending on the scoring system. Normalising to a bit score $S' = (\lambda S - \ln K)/\ln 2$ gives $E = m\,n\,2^{-S'}$ , independent of the scoring matrix. The probability of at least one chance hit is $1-e^{-E}$ . In the demo, $\lambda$ is computed from the scoring rule and $K = 0.46$ is the tabulated value for $+1/-2$ on uniform DNA; on random pairs the predicted $E$ matched our simulation.

Protein search replaces exact words with neighbourhood words: every $w$ -letter word (typically $w=3$ ) whose score against a query word, under a matrix such as BLOSUM62, reaches a threshold $T$ . A higher $T$ is faster and less sensitive, the same trade-off as $w$ for DNA.

The family. Lipman and Pearson's FASTA (1985) used a similar intuition with $k$ -tuples and diagonals. Gapped BLAST and PSI-BLAST (Altschul et al., 1997) required two nearby hits on one diagonal before extending, allowed gaps, and ran about three times faster than the 1990 program. PSI-BLAST iterates, building a position-specific score matrix from the hits to find remote relatives.

Where It Matters

"Find the needle without reading the whole haystack" is the shape of many biological questions:

Function by homology: a new gene that resembles a well-studied one probably does a similar job. This is how newly sequenced genomes get their first annotation.
Identifying organisms: a fragment of DNA from a sample can be matched against reference databases to tell what species, or what pathogen, it came from.
Checking your own data: finding vector or bacterial contamination in a sequencing run.
Remote relatives: iterated searches such as PSI-BLAST reveal protein families whose members share too little sequence for a direct search to notice.

The idea has outgrown its origin. Genome tools such as genome assemblers index sequences by $k$ -mers; later search programs such as DIAMOND and MMseqs2 refine seeding with reduced alphabets and fast vectorised extension; and PatternHunter showed that spaced seeds (patterns with wildcard positions) find more homologs than contiguous words of the same weight. Hash-based sketches (locality-sensitive hashing) go further, giving up alignment entirely in exchange for even faster similarity estimates. The exact method all of them are measured against is the dynamic programming of sequence alignment.

Conclusion

BLAST's insight is that exactness is expensive and most of it is wasted. A short exact word is cheap to find with an index; extending around it with an X-drop rule is cheap to do; and the Karlin–Altschul statistics tell you whether the result could plausibly be a coincidence. What you give up is the guarantee — a distant homolog without a long enough exact word can slip through — and what you gain is a search that finishes while you wait.

The knob you saw in the demo, the word size, is the trade-off every heuristic faces: the more candidates you examine, the more you find and the longer you wait. Decades on, the answer to "what is this sequence?" still starts with seeds, and the exact Smith–Waterman table is what everything else is measured against.

Comments

Loading comments...

https://www.kipuhub.com/en/article/blast-sequence-search/Content licensed under CC BY-NC 4.0.

Cookies & privacy

Seed, extend, and skip almost everything