Sequence Alignment Algorithms — A Historical and Algorithmic Survey

Samuel Seligardi · mentor: Paolo Boldi

The survey traces pairwise sequence alignment from the precursors of the 1960s through Needleman–Wunsch (1970), Hirschberg (1975) and Smith–Waterman (1981) to the seed-and-extend heuristics. This repository holds it together with the code and the measurements behind it; section numbers below refer to the survey.

The survey itself is Survey_Sequence_Alignment_Algorithms.pdf, at the root of this repository.

Three directories:

• algorithms/ — Python implementations of the three dynamic-programming algorithms, one file per variant. Command-line scripts with no dependencies beyond the standard library.
• benchmark/ — the C harness that produced the timings of §7, its Python driver, and the two CSV files those timings live in.
• data/ — the two sequences of the worked example, taken from the 1981 paper.

These are teaching implementations, written to be read next to the paper they come from and next to the pseudocode in the survey.

algorithms/

File	What it does	Paper	Survey
NWv1.py	match reward only — with no gap penalty this is exactly LCS	Needleman & Wunsch (1970)	§4
NWv2.py	adds a linear gap penalty	Needleman & Wunsch (1970)	§4
NWv3.py	adds a mismatch penalty	Needleman & Wunsch (1970)	§4
NWv4.py	match, mismatch and gap all tunable	Needleman & Wunsch (1970)	§4
NWv4_gamma.py	as NWv4, but charging a general gap weight γ(k) per run rather than per column	Needleman & Wunsch (1970)	§6.6
Hirschberg_A.py	the full table — O(mn) time, O(mn) space	Hirschberg (1975)	§5
Hirschberg_B.py	the last row only — O(mn) time, O(m+n) space	Hirschberg (1975)	§5
Hirschberg_C.py	an explicit LCS string by divide-and-conquer — O(mn) time, O(m+n) space	Hirschberg (1975)	§5
SW.py	local alignment with an affine gap weight	Smith & Waterman (1981)	§6

Full references: Needleman & Wunsch, J. Mol. Biol. 48 (1970) 443–453; Hirschberg, Comm. ACM 18 (1975) 341–343; Smith & Waterman, J. Mol. Biol. 147 (1981) 195–197.

The four Needleman–Wunsch versions are incremental by design: each adds one element of the scoring model, so its effect can be seen in isolation. The matrix is filled by anti-diagonal, as the 1970 paper presents it. NWv1.py with gap = 0 is mathematically equivalent to LCS — the score equals the length of the longest common subsequence. That is the link between §4 and §5, and the reason the benchmark scores Needleman–Wunsch as an LCS: the operation count is the same either way.

Hirschberg_A/B/C.py are Hirschberg's own Algorithms A, B and C. B is the building block C calls: C evaluates the middle row in linear space, splits at the optimal crossing point, and recurses.

SW.py implements the 1981 recurrence as written — the two inner maxima rescan every candidate gap length, which is what makes it O(mn(m+n)) rather than O(mn). The zero in the maximum is the whole idea: it discards negative running scores, so an alignment may start anywhere rather than being forced to begin at the first position.

Running them

All nine share one signature, the third argument truncating both inputs (default 100):

python <script> <file1> <file2> [length]

python algorithms/NWv1.py         data/sw81_x.txt data/sw81_y.txt 100
python algorithms/Hirschberg_C.py data/sw81_x.txt data/sw81_y.txt 100

SW.py and NWv4_gamma.py also expose the scoring weights. The defaults are the ones Smith and Waterman use in the paper — match +1, mismatch −1/3, gap 1 + 1/3 per symbol extended:

python algorithms/SW.py data/sw81_x.txt data/sw81_y.txt 100 \
    --match 1 --mismatch -0.33 --gap-base 1 --gap-ext 0.33

Both file arguments are required; omitting them prints the usage message. Run from the repository root, SW.py and NWv4_gamma.py on those two sequences reproduce the local and global alignments compared in §6.6.

data/

Two sequences, and only two. They are the ones Smith and Waterman align in the paper itself (J. Mol. Biol. 147 (1981) 195–197, Figure 1):

sw81_x.txt   AAUGCCAUUGACGG
sw81_y.txt   CAGCCUCGCUUAG

They are here so that the worked example of §6.6 — the local-versus-global comparison — can be regenerated rather than taken on trust:

python algorithms/SW.py         data/sw81_x.txt data/sw81_y.txt   # local
python algorithms/NWv4_gamma.py data/sw81_x.txt data/sw81_y.txt   # global

Under identical scoring, the first yields GCCAUUG against GCC-UCG — seven columns, five matches, 71.4% identity — while the second is forced to place every symbol and yields sixteen columns scoring 1.33, at 50.0% identity. That contrast is the point of §6.6.

Every other input in the survey's benchmark is generated internally, so no further data is needed: alignment_bench.c builds its own sequences from a seeded generator.

benchmark/

alignment_bench.c is a single-threaded C harness that times the dynamic-programming fill of Needleman–Wunsch, Hirschberg and Smith–Waterman over a range of sequence lengths. benchmark.py compiles it, runs it and draws the plots. C rather than the Python above: the constant factor is what puts n = 10⁵ within reach.

Requires gcc on PATH, plus numpy, matplotlib and Pillow.

cd benchmark
python benchmark.py                          # full grid sweep, then the plots
python benchmark.py --diagonal               # m = n only, over a wider range
python benchmark.py --algo nw_lcs sw         # restrict to some algorithms
python benchmark.py --steps heatmap          # redraw from an existing CSV

Phases for --steps are bench heatmap gif in grid mode and bench plot in diagonal mode. The gif phase is opt-in — it draws the timing surface as a rotating 3-D animation, which costs minutes and several MB per algorithm and which no PDF can display. Note that --steps bench writes only the CSV files; the cached .npy grids the plots are drawn from are produced by the heatmap and gif phases. Note also that --algo restricts what the binary prints, and the driver truncates the output file, so re-measuring a single algorithm overwrites the others rather than updating them in place.

To build and run the harness alone:

gcc -O2 -o alignment_bench alignment_bench.c
./alignment_bench > grid_timings.csv

The measurements

grid_timings.csv sweeps m and n independently; diagonal_timings.csv restricts to m = n, where a larger budget reaches n = 10⁵. These two files are what the tables and figures of §7 are built from. An instance whose predicted operation count exceeds a fixed budget is skipped rather than run and recorded as nan.

Measured on an Intel Core i7-13620H (1.25 MB of L2 per performance core, 24 MB of shared L3, 16 GB of RAM) under Windows 11, with gcc 16.1.0 (MinGW-W64) and Python 3.12.4. One run per instance, no repetitions and no warm-up, so small instances carry visible scatter.