DNA Sequencing: Smith-Waterman Alignment and SNP Detection

Sequencing produces fragments, not a genome

Modern short-read sequencers do not read a chromosome start to end. They shred many copies of the DNA into millions of overlapping fragments — "reads" — typically 100 to 300 bases long, and read each one independently and imperfectly, with roughly 0.1 to 1 percent of bases mis-called depending on the platform. Turning that pile of short, error-prone fragments back into a coherent genome sequence, or into a diff against a known reference, is entirely an alignment problem, and alignment is where the actual computer science lives.

Smith-Waterman: finding the best local match

The Smith-Waterman algorithm (1981) is a dynamic-programming method that finds the highest-scoring local alignment between two sequences — the best matching subregion, ignoring the parts that do not correspond. It builds a scoring matrix cell by cell, where each cell represents the best alignment score ending at that pair of positions:

H(i,j) = max(
  0,                                       // reset — this is what makes it LOCAL
  H(i-1, j-1) + score(a[i], b[j]),         // match or mismatch, diagonal
  H(i-1, j)   - gapPenalty,                // deletion, from above
  H(i,   j-1) - gapPenalty                 // insertion, from the left
)

The max(0, …) term is the whole invention: it lets a bad-scoring alignment reset to zero rather than dragging a negative score forward, so the algorithm naturally discovers where a real local match begins. Trace back from the highest-scoring cell in the matrix, following the direction each score came from, until you hit a zero, and you have recovered the best-matching subsequence. Match/mismatch scores typically come from a substitution matrix (a simple +1/-1 scheme for DNA, or something like BLOSUM for proteins), and gaps model insertions or deletions that really do happen in DNA due to replication slippage.

Why not the exact global version?

Smith-Waterman is a one-line modification of the older Needleman-Wunsch algorithm (1970), which finds the best global alignment — forcing the two sequences to align end to end, penalising unmatched overhangs. That is the right tool for comparing two full-length, roughly-the-same-length sequences such as two versions of the same gene. It is the wrong tool for mapping a 150-base read against a 3-billion-base reference genome, where the read corresponds to one small region and everything else is irrelevant — exactly the case Smith-Waterman's local, resettable score handles correctly.

From alignment to variant calling

A single-nucleotide polymorphism (SNP) is a single base that differs from the reference at a given position. Once millions of reads are aligned, calling a SNP is a statistical decision, not a single lookup: a genome sequenced at typical clinical depth is covered by 20 to 30 independent reads at each position, and a real SNP shows the substitution consistently across many of them at a plausible allele fraction (about 50% for one heterozygous copy, close to 100% for a homozygous change). A sequencing error, by contrast, is essentially random per read and rarely reproduces identically across many independent reads of the same base — which is exactly the statistical signal variant callers use to separate real biology from machine noise.

Scaling it up: full-genome alignment in practice

Running full Smith-Waterman — an O(n·m) dynamic program — of a 150-base read against a 3-billion-base genome for every one of hundreds of millions of reads would be computationally impossible. Real aligners (BWA, Bowtie) instead use an FM-index, a compressed full-text index built from the Burrows-Wheeler transform of the reference genome, to instantly narrow a read down to a handful of candidate locations in roughly O(read length) time. Full Smith-Waterman-style scoring is then only run locally, around those few candidate positions, to get a precise, gapped alignment and quality score — combining a near-exact-match index for speed with dynamic programming for accuracy exactly where it is needed.

Frequently asked questions

Try it live

Everything above runs in your browser — open DNA Sequencing and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.