How Read Aligners Work: Mapping Reads to a Reference Genome

Every variant call, every expression count, every peak in a ChIP-seq track starts with the same unglamorous question: where on the genome did this read come from? Answering that question a few hundred million times per sample, in minutes rather than days, is what read aligners like BWA-MEM, Bowtie2, and HISAT2 are built to do. This post opens up the machinery underneath the bwa mem or bowtie2 command you run every day.

Why We Need an Index At All

The naive way to find a read's origin is brute force: slide the read along the reference and score the match at every position, the way Smith-Waterman local alignment works. For one 150 bp read against a small reference that is fine. For a human genome of roughly 3.2 billion bases and a lane of 400 million reads, brute force is computationally out of reach — each read would cost work proportional to the length of the genome, repeated hundreds of millions of times.

The fix is to do the expensive work exactly once. Before any reads are aligned, the aligner builds an index of the reference genome — a structure that lets it ask "where does this sequence occur?" in time proportional to the length of the query, not the genome. That is why every aligner ships a separate index command (bwa index, bowtie2-build, hisat2-build) run once per reference and reused for every sample afterward.

The Suffix Array Intuition

The classic way to build a searchable index of a string is a suffix array: list every suffix of the reference (the substring starting at each position and running to the end), then sort those suffixes alphabetically. Once sorted, every occurrence of a query pattern in the genome corresponds to a contiguous block in that sorted list, because all suffixes that start with the same prefix sit next to each other after sorting. Finding a query then becomes a binary search over the sorted suffixes — O(m log n) for a query of length m against a reference of length n, instead of scanning the whole genome.

The catch is memory. A full suffix array stores one integer position per base of the genome, so for a 3+ billion base human genome the array alone runs several times larger than the raw sequence itself. Something more compact was needed before genome-scale short-read alignment could actually run on ordinary hardware — that something is the Burrows-Wheeler Transform.

The Burrows-Wheeler Transform: Compressing the Search

The Burrows-Wheeler Transform, or BWT, is a reversible rearrangement of a string that was originally invented for data compression, decades before anyone used it for sequence alignment. Conceptually, you form every rotation of the reference sequence, sort those rotations alphabetically (this sort is closely related to the suffix array sort above), and then read off the last character of each sorted rotation. That final column is the BWT.

Two properties make it valuable here. First, because the rotations are sorted, the BWT tends to group identical or highly repetitive contexts together, which is exactly the kind of run-length structure that compresses well — the same property bgzip-compressed BAM files exploit downstream. Second, and more importantly for alignment, the BWT is derived from the same sorted order as the suffix array, so it inherits the ability to answer "does this pattern occur, and where" — without needing to store the array itself.

The FM-Index: Searching Without Decompressing

A bare BWT string alone cannot be searched quickly; you would still need to decompress or scan it. The FM-index, introduced by Ferragina and Manzini, pairs the BWT with a small amount of auxiliary bookkeeping — checkpointed counts of how many times each base has appeared up to a given position — that together support backward search. Starting from the last character of the query and working toward the first, the FM-index narrows a range in the sorted-rotation order one character at a time, each step costing roughly constant time regardless of genome size.

This is the innovation that made Bowtie and BWA practical. Heng Li and Richard Durbin's 2009 paper on BWA describes exactly this approach, and it is why a full index of the human genome fits in only a few gigabytes of RAM instead of the tens of gigabytes a raw suffix array would need — suffix-array-quality search at a fraction of the memory, small enough for a laptop.

Seed-and-Extend: From Exact Matches to Real Alignments

Backward search over an FM-index is powerful, but it only finds exact matches. Real sequencing reads are not exact copies of the reference — they carry sequencing errors, real SNPs, and small insertions or deletions. If you demanded a perfect match for the whole read, most reads with even one mismatch would fail to align at all.

The solution used by every mainstream aligner is seed-and-extend, and it has two stages:

1. Seed. Break the read into short exact-match anchors and locate each one in the genome using the FM-index. BWA-MEM uses maximal exact matches (MEMs) — the longest exact substrings of the read that occur in the reference — while Bowtie2 and HISAT2 extract multiple shorter seeds from different offsets along the read. Seeding is fast precisely because it only needs exact matching.
2. Extend. Each seed's genomic location is a candidate anchor. Around every anchor, the aligner runs a slower but far more forgiving dynamic-programming alignment — a Smith-Waterman-style algorithm with affine gap penalties — to check whether the seed actually grows into a full, high-scoring alignment that tolerates mismatches and small indels. Only the handful of candidate regions near real seed hits get this expensive treatment, not the whole genome.

That division of labor is the whole trick: the FM-index cheaply discards almost the entire genome as irrelevant, and dynamic programming — accurate but slow — only ever runs on the handful of places that already look promising.

BWA-MEM, Bowtie2, and HISAT2: Same Ideas, Different Jobs

All three tools sit on the BWT/FM-index and seed-and-extend foundation, but they are tuned for different data.

The splice-aware column is the key distinction for RNA-seq work. A read spanning an exon-exon junction has no contiguous match anywhere in the genome, because the intron between the exons is absent from the mature transcript. HISAT2 solves this with a hierarchical graph FM-index: alongside the whole-genome index, it layers many small local indexes around splice sites so a read can be seeded and extended across a junction it would otherwise miss entirely. This is also why splice-aware alignment belongs upstream of the kind of differential expression workflow built on featureCounts and DESeq2 — a genomic (non-splice-aware) aligner like BWA-MEM will systematically fail to place exon-spanning reads, undercounting genes with short exons.

Mapping Quality (MAPQ): How Confident Is This Placement?

Finding an alignment is not the same as finding the right one. Repetitive elements, gene duplications, and short reads all create places where a read aligns equally well to two or more genomic locations. To flag that ambiguity, every aligner reports a mapping quality score, MAPQ, in column 5 of the SAM/BAM record.

MAPQ is a Phred-scaled estimate of the probability that the reported position is wrong, so higher is better, capped at an aligner-specific maximum — commonly 60 for BWA-MEM. It is computed largely from the gap between the best-scoring alignment and the next-best one found during the extend step: one clearly dominant alignment gets a high MAPQ, while a read scoring nearly identically at two or more locations gets a low MAPQ, often 0 or 1, regardless of how good the alignment score itself looks.

This matters downstream. Variant callers typically discard reads below a MAPQ threshold — filtering to MAPQ >= 20 or MAPQ >= 30 is common — because a confidently wrong placement is worse than no placement at all for calling a SNP or indel. For how those filtered alignments become variant calls, see our guide to GATK4.

Practical Takeaways

• Indexing exists to move the expensive work out of the per-read loop: build the FM-index once with bwa index or hisat2-build, then align millions of reads against it cheaply.
• The FM-index gives you suffix-array-quality exact search in a fraction of the memory, which is why genome-scale alignment fits on a laptop rather than requiring a cluster.
• Seed-and-extend is the pattern to remember: fast exact seeding via the FM-index narrows the genome down to a handful of candidates, and slow, accurate dynamic programming only runs on those candidates.
• Choose the aligner for the biology, not just the algorithm: BWA-MEM and Bowtie2 assume contiguous genomic matches, while HISAT2's graph FM-index is built specifically to span exon-exon junctions in RNA-seq data.
• Always check MAPQ, not just whether a read aligned. A high alignment score in a repetitive region can still carry a MAPQ near zero, and downstream tools should be filtering on it.

Understanding this pipeline also explains why installing these tools correctly matters — if you have not set one up yet, our guide to conda environments and Bioconda channels covers the isolated-environment approach used across every aligner install on this site.