Genomics

The central dogma. The flow of genetic information in a cell follows a single direction: DNA is transcribed into RNA, which is translated into protein. DNA is a double-stranded polymer made of four nucleotides - adenine (A), cytosine (C), guanine (G), thymine (T) - held together by complementary base pairing (A–T, C–G). The complete DNA sequence of an organism is its genome. Genomics is the study of genomes: their structure, content, and evolution.

Computational challenges. A human genome contains roughly $3 \times 10^9$ base pairs. No current technology reads it in one pass: sequencers produce millions of short reads (fragments of 100–30,000 bp depending on the platform), and computational methods must reconstruct the original sequence, identify variants, and annotate functional elements. The core algorithmic problems are sequence alignment (comparing two sequences) and genome assembly (reconstructing a genome from overlapping reads).

String algorithms at the core. A DNA sequence is a string over the alphabet $\Sigma = \{\texttt{A}, \texttt{C}, \texttt{G}, \texttt{T}\}$. Genomics imports and adapts classical string algorithms - dynamic programming, suffix arrays, hash tables - and scales them to gigabase inputs. The dominant complexity constraints are time (reads number in the millions) and memory (a genome-scale suffix array requires tens of gigabytes).

The sequencing problem. A sequencer cannot read a whole chromosome in one pass. Instead, it produces a large number of short reads: substrings sampled (approximately at random) from the genome. Reconstruction of the original genome from these reads is the assembly problem. The key parameters of a sequencing technology are: read length, throughput (total bases per run), error rate, and error profile (substitutions vs. insertions/deletions).

Short-read sequencing (Illumina). The dominant technology since $\sim$2007. DNA is fragmented, adapters are ligated to both ends, and fragments are amplified on a flow cell by bridge PCR, forming clusters. Sequencing proceeds by synthesis: at each cycle, fluorescently labeled reversible-terminator nucleotides are incorporated, imaged, and cleaved. This yields reads of 100–300 bp with very low error rates ($\sim$0.1%), primarily substitution errors. A single run produces hundreds of gigabases of data. The short read length is the main limitation: it makes it hard to assemble repetitive regions longer than a read.

Long-read sequencing (PacBio, Oxford Nanopore). Two platforms produce reads of 10,000–1,000,000 bp, enabling assembly across repetitive regions. PacBio (SMRT sequencing) observes polymerase activity in zero-mode waveguides; Oxford Nanopore measures ionic current changes as DNA passes through a protein pore. Both have higher error rates ($\sim$1–15%) concentrated in homopolymer runs (consecutive identical bases), but sequencing errors are largely random and can be corrected by high coverage or hybrid assembly with short reads.

Sequence file formats: FASTA and FASTQ. Sequencing data and reference genomes are exchanged in two plain-text formats. FASTA stores sequences without quality information: each entry is a header line beginning with > followed by an arbitrary identifier, then the sequence on subsequent lines. FASTQ stores reads with per-base quality scores: four lines per read - a header beginning with @, the sequence, a separator line +, and one ASCII-encoded quality character per base. The Phred quality score at a base is $Q = -10 \log_{10} P_{\text{err}}$, where $P_{\text{err}}$ is the estimated probability that the base call is incorrect, encoded in the file as the character $\texttt{chr}(Q + 33)$. So $Q = 20$ corresponds to error probability $10^{-2}$, $Q = 30$ to $10^{-3}$, and $Q = 40$ to $10^{-4}$. Quality scores are essential for filtering low-confidence bases and for weighting evidence in variant calling.

Coverage intuition. Think of tearing a book into short, randomly placed strips. To reconstruct the text you need enough strips so that every word is covered by at least one strip. With only one strip per region you cannot distinguish a smudged word (sequencing error) from the real text; with 30 overlapping strips showing the same word you can be confident. Two complications arise in practice. First, coverage is uneven: some regions are over-represented (biased amplification, high GC content) and others are missed entirely. Second, if the same sentence is repeated in many chapters, strips from those regions are indistinguishable and you cannot tell where each one belongs. These two phenomena -- non-uniform coverage and repeats -- are why the idealized Poisson model below is only an approximation, and why 30$\times$ coverage is needed in practice rather than the $\approx 5\times$ that the model predicts.

The pattern matching problem. The simplest string query in genomics: given a short pattern $p$ of length $m$ and a long text $t$ of length $n$, find all positions $i$ such that $t[i..i+m-1] = p$. Practical examples include locating restriction enzyme cut sites ($m \approx 6$), checking whether a primer binds to a reference ($m \approx 20$), and finding exact seed matches as a subroutine of read alignment. Because reads may originate from either strand of the double helix, one must always search for both $p$ and its reverse complement $\overline{p}$.

Naïve matching. The naïve algorithm slides a window of length $m$ over $t$ and compares all characters at each position. It runs in $O(nm)$ time and $O(1)$ extra space. For a human genome ($n = 3 \times 10^9$, $m = 20$), this is $6 \times 10^{10}$ comparisons. Early termination on mismatch makes the average case much faster: for a uniformly random DNA string, the expected number of comparisons per window is $O(1)$ since the first character mismatches with probability $3/4$. In the worst case ($t = \texttt{AAA}\cdots$, $p = \texttt{AAA}\cdots\texttt{B}$) every window requires $m$ comparisons. This worst case motivates the KMP and Boyer-Moore algorithms, which guarantee $O(n)$ time.

def rc(seq):
    return seq.translate(str.maketrans('ACGT', 'TGCA'))[::-1]

def exact_match(pattern, text):
    m, n = len(pattern), len(text)
    return [i for i in range(n - m + 1) if text[i:i+m] == pattern]

def match_both_strands(pattern, text):
    return exact_match(pattern, text), exact_match(rc(pattern), text)

Indexed search: BWT and FM-index. Scanning the full text for every query is wasteful when the text is fixed (a reference genome) and many queries (millions of reads) will be issued against it. The idea of an index is to preprocess the text once so that subsequent queries take time depending only on the pattern length $m$, not on the text length $n$.

The Burrows–Wheeler transform (BWT) of a string $t$ (terminated by a sentinel symbol $\#$ smaller than every alphabet character) is the last column of the matrix whose rows are all cyclic rotations of $t\#$ sorted lexicographically. Remarkably, $t$ can be recovered exactly from its BWT alone, and together with two small auxiliary tables -- the $C$ array (number of characters in $t$ smaller than each symbol) and rank counts -- the BWT supports backward search: given a pattern $p$ of length $m$, the range of suffix-array positions matching $p$ is computed in $O(m)$ time, independent of $n$. This is the FM-index (Ferragina–Manzini, 2000).

Modern short-read aligners (BWA, Bowtie) build the FM-index of the reference once ($\sim$3 GB for the human genome) and align each read with a handful of BWT operations. The FM-index combines the query speed of a suffix tree with a memory footprint close to that of the raw text.

def hamming(a, b):
    assert len(a) == len(b)
    return sum(x != y for x, y in zip(a, b))

Why dynamic programming? To align $a$ (length $m$) against $b$ (length $n$) we need to decide, for each character, whether it is matched to a character of the other string or to a gap. The number of possible alignments is exponential: roughly $\binom{m+n}{m}$, which grows faster than $2^m$ for comparable lengths. Enumerating them all is hopeless.

The saving observation is optimal substructure: the last column of any alignment of $a[1..i]$ with $b[1..j]$ is one of exactly three things.

1. $a_i$ is aligned with $b_j$ (a match or mismatch).
2. $a_i$ is aligned with a gap on the $b$ side.
3. $b_j$ is aligned with a gap on the $a$ side.

In every case, removing that last column leaves an alignment of a strictly shorter prefix pair: $(a[1..i-1], b[1..j-1])$, $(a[1..i-1], b[1..j])$, or $(a[1..i], b[1..j-1])$ respectively. If the original alignment is optimal, the leftover prefix alignment must also be optimal for its own prefix pair (otherwise we could swap in a better prefix alignment and improve the total score -- a contradiction). This is the exchange argument that justifies recursion.

So define $F(i, j)$ as the score of the best alignment of $a[1..i]$ with $b[1..j]$. We can compute $F(i,j)$ from three smaller subproblems, building the table from top-left to bottom-right. Computing the whole $(m+1)\times(n+1)$ table takes $O(mn)$ time and space.

def needleman_wunsch(a, b, match=1, mismatch=-1, gap=-1):
    m, n = len(a), len(b)
    dp = [[i * gap] + [0] * n for i in range(m + 1)]
    dp[0] = [j * gap for j in range(n + 1)]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            sub = match if a[i-1] == b[j-1] else mismatch
            dp[i][j] = max(
                dp[i-1][j-1] + sub, dp[i-1][j] + gap, dp[i][j-1] + gap
            )
    ra, rb, i, j = [], [], m, n
    while i or j:
        sub = match if i and j and a[i-1] == b[j-1] else mismatch
        if i and j and dp[i][j] == dp[i-1][j-1] + sub:
            ra.append(a[i-1])
            rb.append(b[j-1])
            i -= 1
            j -= 1
        elif i and dp[i][j] == dp[i-1][j] + gap:
            ra.append(a[i-1])
            rb.append('-')
            i -= 1
        else:
            ra.append('-')
            rb.append(b[j-1])
            j -= 1
    return ''.join(reversed(ra)), ''.join(reversed(rb)), dp[m][n]

From global to local alignment. Needleman–Wunsch aligns all of $a$ against all of $b$. End gaps are penalized just like interior gaps. This is the right model when two sequences are known to be related along their full length (e.g. two homologous genes of similar size).

For many biological questions the situation is different: $a$ might be a short query (a read, a protein domain) and $b$ a long database sequence, or both sequences may share only a conserved region flanked by unrelated sequence. In these cases, forcing a full alignment penalizes the unrelated flanking regions heavily and yields a meaningless score. What we want is to find the pair of substrings $(a[i..i'],\, b[j..j'])$ whose alignment score is highest -- ignoring everything outside them.

The trick to turning NW into a local aligner is simple: we allow the alignment score to "restart" at any position by adding $0$ as a fourth option in the recurrence. Concretely, define $H(i,j)$ as the score of the best alignment ending at $a_i$ and $b_j$, or $0$ if no alignment ending here has a positive score. Then extending from a cell with $H = 0$ is the same as starting a fresh alignment at that point. The best local alignment score is $\max_{i,j} H(i,j)$, found anywhere in the table, not just at the bottom-right corner. Traceback starts at that maximum and stops as soon as it reaches a cell with $H = 0$ -- that cell marks the start of the local alignment.

def smith_waterman(a, b, match=2, mismatch=-1, gap=-1):
    m, n = len(a), len(b)
    dp = [[0] * (n + 1) for _ in range(m + 1)]
    best, pos = 0, (0, 0)
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            sub = match if a[i-1] == b[j-1] else mismatch
            dp[i][j] = max(
                0, dp[i-1][j-1] + sub, dp[i-1][j] + gap, dp[i][j-1] + gap
            )
            if dp[i][j] > best:
                best, pos = dp[i][j], (i, j)
    ra, rb = [], []
    i, j = pos
    while dp[i][j]:
        sub = match if a[i-1] == b[j-1] else mismatch
        if dp[i][j] == dp[i-1][j-1] + sub:
            ra.append(a[i-1])
            rb.append(b[j-1])
            i -= 1
            j -= 1
        elif dp[i][j] == dp[i-1][j] + gap:
            ra.append(a[i-1])
            rb.append('-')
            i -= 1
        else:
            ra.append('-')
            rb.append(b[j-1])
            j -= 1
    return ''.join(reversed(ra)), ''.join(reversed(rb)), best

Affine gap penalties. The linear gap penalty $-d$ per gap column (used above) treats a gap of length $\ell$ as $\ell$ independent gaps. Biologically this is wrong: a single insertion of 10 bases is one mutational event, far more likely than 10 separate single-base insertions. An affine gap penalty charges a fixed cost $-o$ to open a gap plus a smaller cost $-e$ to extend it by one additional base, so a gap of length $\ell$ costs $o + (\ell - 1) e$ with $o \gg e$ (typical values: $o = 10$, $e = 1$).

Implementing affine penalties requires tracking whether the previous column was already inside a gap, which the single-table NW/SW recurrence cannot do. The fix is three tables: $M(i,j)$ for alignments ending in a match/mismatch, $X(i,j)$ for alignments ending in a gap in $b$ (consuming a character of $a$), and $Y(i,j)$ for alignments ending in a gap in $a$. Each table's recurrence distinguishes opening a gap (transition from $M$, costing $-o$) from extending a gap (staying in $X$ or $Y$, costing $-e$). The total complexity remains $O(mn)$ in both time and space. This is the Gotoh algorithm (1982); it is the default in modern aligners (BWA, minimap2) and BLAST.

From exact to heuristic alignment: BLAST. Needleman–Wunsch and Smith–Waterman are $O(mn)$. When $m = 300$ (a read) and $n = 3 \times 10^9$ (the human genome), this is $9 \times 10^{11}$ operations per read - far too slow. BLAST (Basic Local Alignment Search Tool) accelerates local alignment by a seed-and-extend heuristic:
(1) find all exact matches of length $W$ (seeds) between the query and the database,
(2) extend each seed in both directions as long as the local alignment score stays above a threshold. This misses alignments with no exact $W$-mer match, but in practice very few biologically meaningful alignments are missed when $W = 11$. Modern short-read aligners (BWA, Bowtie2) use different data structures (FM-index, Burrows–Wheeler transform) to achieve $O(n + m)$ alignment.

The assembly problem. Given a set of reads $R = \{r_1, \ldots, r_N\}$ sampled from an unknown genome $G$, reconstruct $G$. This is fundamentally ambiguous: repetitive regions of length greater than a read cannot be resolved from the reads alone. In practice, the output is not a single sequence but a set of contigs (contiguous assembled sequences) and scaffolds (contigs linked by mate-pair information, with gaps of estimated size). Assembly quality is measured by the N50: the length $L$ such that contigs of length $\geq L$ cover at least half the total assembled sequence.

def overlap(r1, r2, min_len):
    return next((k for k in range(min(len(r1), len(r2)), min_len - 1, -1)
                 if r1[-k:] == r2[:k]), 0)

def overlap_graph(reads, min_len):
    return {
        i: [(j, k) for j, r2 in enumerate(reads)
            if i != j and (k := overlap(r1, r2, min_len)) > 0]
        for i, r1 in enumerate(reads)
    }

The de Bruijn idea: switching from reads to k-mers. The overlap graph approach has two problems. First, computing all pairwise overlaps between $N$ reads costs $O(N^2)$, which is prohibitive for millions of reads. Second, finding a Hamiltonian path is NP-hard.

The de Bruijn graph sidesteps both by changing what the vertices and edges represent. Instead of tracking reads, we track k-mers. The key observation is: two consecutive k-mers in a sequence necessarily share a $(k-1)$-mer overlap -- the suffix of the first equals the prefix of the second. For example, with $k = 3$: the k-mers ACG and CGT share the overlap CG.

We exploit this as follows: make each distinct $(k-1)$-mer a node, and represent each k-mer as a directed edge from its length-$(k-1)$ prefix to its length-$(k-1)$ suffix. The genome then corresponds to a path that traverses every edge (k-mer) exactly once -- an Eulerian path. Finding an Eulerian path takes $O(|E|)$ time (Hierholzer's algorithm), a dramatic improvement over the NP-hard Hamiltonian path problem.

Recovering the sequence from the path is immediate: consecutive nodes in the path share a $(k-2)$-mer overlap, so the assembled sequence is the first node concatenated with the last character of each subsequent node (see the assemble function below).

from collections import Counter, defaultdict

def debruijn_graph(reads, k):
    graph = defaultdict(list)
    for read in reads:
        for i in range(len(read) - k + 1):
            graph[read[i:i+k-1]].append(read[i+1:i+k])
    return graph

def eulerian_path(graph):
    in_deg = Counter(v for vs in graph.values() for v in vs)
    out_deg = {u: len(vs) for u, vs in graph.items()}
    start = next(
        (u for u in graph if out_deg[u] - in_deg[u] == 1), next(iter(graph))
    )
    adj = {u: list(vs) for u, vs in graph.items()}
    stack, path = [start], []
    while stack:
        if adj.get(stack[-1]):
            stack.append(adj[stack[-1]].pop())
        else:
            path.append(stack.pop())
    return path[::-1]

def assemble(reads, k):
    path = eulerian_path(debruijn_graph(reads, k))
    return path[0] + ''.join(node[-1] for node in path[1:])

What is annotation. A raw genome sequence is just a string. Genome annotation is the process of identifying the location and function of its elements: protein-coding genes, non-coding RNAs, regulatory sequences, and repetitive elements. It proceeds in two phases: structural annotation (where are the genes?) and functional annotation (what do they do?). The two main strategies for gene finding are ab initio prediction (from sequence signals alone) and homology-based prediction (by similarity to known genes).

The genetic code and reading frames. The information in a gene is read in consecutive non-overlapping triplets called codons. Each codon specifies one amino acid in the resulting protein, or signals where to start and stop translation. There are $4^3 = 64$ possible codons but only 20 amino acids; most amino acids are encoded by multiple codons (the code is redundant). Two special roles:

• Start codon: ATG (encodes methionine). Translation begins here.
• Stop codons: TAA, TAG, TGA. They encode no amino acid; they terminate translation.

Since codons are triplets, there are three possible ways to partition a sequence into consecutive triplets depending on where you start: positions $\{0, 3, 6, \ldots\}$, $\{1, 4, 7, \ldots\}$, or $\{2, 5, 8, \ldots\}$. Each choice is called a reading frame. Since a gene can lie on either strand of the double helix, there are $3 \times 2 = 6$ reading frames to consider for any given sequence.

STOP_CODONS = {'TAA', 'TAG', 'TGA'}

def find_orfs(seq, min_codons=100):
    orfs = []
    for frame in range(3):
        codons = [seq[i:i+3] for i in range(frame, len(seq) - 2, 3)]
        starts = [i for i, c in enumerate(codons) if c == 'ATG']
        for start in starts:
            stop = next(
                (i for i in range(start + 1, len(codons))
                 if codons[i] in STOP_CODONS),
                None,
            )
            if stop and stop - start >= min_codons:
                orfs.append((frame + start * 3, frame + (stop + 1) * 3, frame))
    return orfs

def reverse_complement(seq):
    return seq.translate(str.maketrans('ACGT', 'TGCA'))[::-1]

def find_all_orfs(seq, min_codons=100):
    rc = reverse_complement(seq)
    reverse = [(len(seq) - e, len(seq) - s, frame + 3)
               for s, e, frame in find_orfs(rc, min_codons)]
    return find_orfs(seq, min_codons) + reverse

Ab initio gene prediction. In eukaryotes, gene structure is complex: exons (coding), introns (non-coding, spliced out), UTRs, and splice signals. Ab initio predictors (Augustus, GeneMark) model these signals using hidden Markov models (HMMs). States represent genomic features (exon, intron, intergenic, splice donor, splice acceptor), and emission probabilities capture codon usage bias and splice consensus sequences. The Viterbi algorithm finds the most likely state sequence - and hence the most likely gene structure - in $O(n \cdot |S|^2)$ time, where $|S|$ is the number of states.

Functional annotation. Once gene locations are predicted, each gene product is assigned a function by similarity search against curated databases: UniProt/Swiss-Prot (manually curated protein sequences), Pfam (protein domain families modeled as profile HMMs), and the Gene Ontology (GO) - a controlled vocabulary of biological process, molecular function, and cellular component terms. A BLAST or HMMer search transfers annotations from known homologs to the new gene. The reliability of transferred annotations decreases as sequence identity drops below $\sim$30%.

The variant calling pipeline. The standard pipeline has three stages.
(1) Read mapping: align reads to the reference genome using a short-read aligner (BWA-MEM, Bowtie2). Each read is placed at its best matching position; the result is a BAM file - a sorted, indexed binary alignment.
(2) Pileup: for each reference position, collect all reads covering that position and count the observed bases.
(3) Genotyping: apply a statistical model to decide whether an observed deviation from the reference is a true variant or a sequencing error. The industry-standard tool is GATK HaplotypeCaller, which uses a local assembly step to handle indels cleanly.

Diploid genomes and alleles. Humans (and most animals) are diploid: every cell carries two copies of each chromosome, one inherited from each parent. At any given genomic position, the two copies may carry the same base (the individual is homozygous) or different bases (the individual is heterozygous). Call the base in the reference genome $R$ and any differing base $A$ (alternative). A genotype at a position is one of $RR$, $RA$, or $AA$. Each sequenced read comes from one of the two chromosomal copies chosen at random with equal probability.

from math import log

def genotype_likelihoods(n, k, error=1e-3):
    return {
        'RR': (n-k) * log(1 - error) + k * log(error),
        'RA': n * log(0.5),
        'AA': k * log(1 - error) + (n-k) * log(error),
    }

def call_genotype(n, k, error=1e-3):
    lls = genotype_likelihoods(n, k, error)
    return max(lls, key=lls.get)

Goals of comparative genomics. Comparing genomes across species reveals which regions are evolutionarily conserved (and therefore likely functional), infers the evolutionary history of genes and species, and identifies gene gains, losses, and rearrangements. The two central problems are multiple sequence alignment (aligning more than two sequences simultaneously) and phylogenetic tree inference (reconstructing the evolutionary tree from sequence data).

Why raw sequence divergence underestimates evolutionary distance. Suppose we compare two sequences and find that 30% of sites differ. The naive estimate of evolutionary distance is "0.3 substitutions per site." But this is too low, for two reasons.

First, back mutations: a site might mutate from A to T and then back to A. Those two mutations are invisible to us because the site looks identical, yet two events occurred. Second, multiple hits: a site mutates from A to T, then from T to G. We observe a difference (A vs G) and count it as one event, but two actually occurred. Both effects cause the observed fraction of differing sites to saturate as sequences diverge further -- the curve flattens toward the random maximum of 75% (since at equilibrium with equal frequencies, 4 equally probable bases, 3 out of 4 random pairs will differ). The true number of mutations is unbounded.

A correction model treats each site as a Markov chain over the four bases and derives the expected fraction of identical sites as a function of the true number of mutations per site $d$. Inverting that function gives a corrected distance estimate from the observed divergence.

from math import log, inf

def jukes_cantor(s1, s2):
    p = sum(a != b for a, b in zip(s1, s2)) / len(s1)
    return -0.75 * log(1 - 4 * p / 3) if p < 0.75 else inf

def upgma(seqs):
    n = len(seqs)
    clusters = {i: [i] for i in range(n)}
    dist = {
        (i, j): jukes_cantor(seqs[i], seqs[j])
        for i in range(n)
        for j in range(i+1, n)
    }
    merges, next_id = [], n
    while len(clusters) > 1:
        i, j = min(dist, key=dist.get)
        merges.append((i, j, dist[i, j] / 2, next_id))
        new = clusters[i] + clusters[j]
        for k in set(clusters) - {i, j}:
            dist[min(next_id, k), max(next_id, k)] = (
                len(clusters[i]) * dist[min(i, k), max(i, k)] +
                len(clusters[j]) * dist[min(j, k), max(j, k)]
            ) / len(new)
        dist = {p: d for p, d in dist.items() if i not in p and j not in p}
        del clusters[i], clusters[j]
        clusters[next_id] = new
        next_id += 1
    return merges