6 Read Alignment: Concepts

6.1 The problem

A sequencing run hands you hundreds of millions of short strings and no idea where any of them came from. Read alignment is the problem of deciding, for each of those reads, which position in a reference genome it was most likely sequenced from — and how it differs from the reference there.

Almost everything downstream depends on getting this right. A variant caller can only see substitutions that survived alignment. An expression estimate is a count of reads assigned to features. Alignment is where the data stops being strings and starts being coordinates, and errors introduced here are very hard to detect later, because everything after it treats the coordinates as given.

Two words that get used interchangeably and should not be. Mapping is finding where a read belongs. Alignment is describing how it matches once it gets there — the exact correspondence of bases, including mismatches and gaps. Most tools do both, which is why the distinction blurs, but they can fail independently: a read can be mapped to the right locus with a bad base-level alignment, and a read with an internally consistent alignment can be at the wrong locus entirely.

6.2 Three questions, three kinds of alignment

Before any algorithm, you have to know which question you are asking.

Global alignment forces both sequences to align end to end. This is the Needleman–Wunsch formulation (Needleman and Wunsch 1970), and it is the right question for two sequences you believe are homologous over their whole length — two versions of the same gene, say.

Local alignment finds the best-matching subsequence of each, ignoring the rest. This is Smith–Waterman (Smith and Waterman 1981), and it is what you want when a short read may only partly correspond to the reference: adapter still attached, a structural breakpoint mid-read, the end of a chromosome.

Semi-global alignment — sometimes “glocal” — requires one sequence to align end to end while letting the other extend past it freely. This is the one short read alignment actually wants most of the time: the read should be consumed entirely, but the genome obviously continues in both directions.

Choosing wrong is a quiet failure. Running an end-to-end alignment on reads that still carry adapter forces the aligner to explain sequence that has no genomic origin, and it will do so, by inventing mismatches near the read’s end.

6.3 Why the exact answer does not scale

All three are solved exactly by dynamic programming, and the cost of that is proportional to the product of the two lengths. For one read against one gene, this is nothing. For a hundred million reads against three billion bases, it is hopeless by many orders of magnitude — and no constant-factor engineering closes a gap that size.

So every practical aligner makes the same structural bet, called seed and extend: find short stretches that match exactly, or nearly so, and use them to nominate a small number of candidate positions. Then run the expensive exact algorithm only in those windows.

The bet is that a true alignment contains at least one exact short match. It is a good bet — a 100 bp read at 1% divergence will contain long exact stretches — and it is the source of essentially every failure mode in the field. Reads that are too divergent, too short, or too repetitive break the assumption, and the aligner does not fail loudly. It reports something.

6.4 The two index families

Seeding needs an index that answers “where does this short string occur?” quickly. Two families dominate, and the choice between them explains most of the differences between the tools in the next chapter.

Suffix-based indexes — the Burrows–Wheeler transform and the FM-index (Ferragina and Manzini 2000) — support exact substring search in time proportional to the query length, independent of how large the reference is, in space close to the compressed reference. Applying this to read alignment (Li and Durbin 2009) is what made whole-genome short-read alignment routine. Their strength is exact, complete answers: every occurrence, guaranteed. Their weakness is that mismatches are expensive, since each one forces the search to branch.

Sketch-based indexes store only a sample of the reference’s k-mers. The standard sampling scheme is the minimizer (Roberts et al. 2004): in each window of consecutive k-mers, keep only the one with the smallest hash value. Two sequences that overlap will tend to select the same k-mer in the shared region, so the sample stays informative while shrinking the index substantially. Sketch indexes tolerate divergence far better and are much faster to query, at the cost of being approximate — a real match can be missed if its minimizers happen not to be sampled.

Why long reads changed the answer

At 1% error, an exact-match index is ideal: exact seeds are long and plentiful. At 10% error, exact seeds of useful length essentially stop existing, and an FM-index spends all its time branching. This is why the long-read era did not just re-tune the short-read aligners — it moved the field to sketch-based seeding plus chaining, where many short, individually unconvincing anchors are assembled into a colinear run whose collective evidence is strong. The concept chapter’s point is that this was a change of question, not a change of implementation.

6.5 Mapping quality is not alignment score

This is the single most misread field in the output, so it gets its own section.

An alignment score says how well the read matches at the position it was placed. A mapping quality (MAPQ) says how confident the aligner is that the position is the right one. It is Phred-scaled:

\[ \text{MAPQ} = -10 \log_{10} \Pr\{\text{the mapping position is wrong}\} \]

The two come apart in both directions, and both directions matter.

A read can align perfectly — a flawless score — and carry MAPQ 0, because it aligns equally perfectly to four other places. A read in a segmental duplication is the standard case. Conversely a read can align poorly, with several mismatches, and still carry a high MAPQ, because nowhere else in the genome is remotely as good.

Watch out

Filtering on MAPQ discards reads from repetitive regions by construction. That is often what you want, and it is never neutral: it silently makes paralogous gene families, immune loci and centromeric sequence disappear from your analysis. If a downstream result concerns a repetitive region, MAPQ filtering is part of your method and belongs in the write-up.

MAPQ is also not comparable across tools. It is a model-dependent estimate, and different aligners scale it differently; a threshold tuned for one is not meaningful for another.

6.6 The output is a contract, not a file format

Alignments are written as SAM, or its compressed binary form BAM, or the reference-relative CRAM (Li et al. 2009). Each record carries the read, the position it was assigned, its MAPQ, a bitwise flag describing pairing and orientation, and a CIGAR string spelling out the alignment as a run-length sequence of matches, insertions, deletions, clips and skips.

The part worth internalizing is that the file’s type is more than its extension. A BAM may or may not be coordinate-sorted; it may or may not have an index alongside it; its coordinates are meaningless without knowing which reference build produced them. Nearly every downstream tool requires a coordinate-sorted, indexed BAM against a known reference, and a plain BAM does not satisfy that requirement.

Think of “sorted, indexed BAM against build X” as the actual type of the data, and “BAM” as merely its encoding. Most pipeline failures in practice are type errors in this sense — a step handed a file that has the right extension and the wrong properties. The failure is rarely a crash; more often it is an empty result or a silently partial one.

Two CIGAR operations are worth separating now, because they are routinely confused. Soft clipping (S) means the aligner did not use the read’s end but kept the sequence in the record. Hard clipping (H) means it discarded it. And N, the skip operation, is not a deletion — it means the read spans a region absent from the transcript, which is how introns are represented.

6.7 Spliced alignment is a different problem

Which brings us to the case that breaks the framework above. A read from a mature transcript may cross an exon boundary, so its two halves align to positions separated by thousands of bases of intron that are genuinely not in the molecule.

An ordinary aligner can only describe that as a deletion, and its gap penalties make a deletion of that size effectively impossible. The read gets clipped, or placed wrongly, or dropped. Spliced aligners treat these as a distinct operation with its own cost model, informed by the fact that introns are not arbitrary: they overwhelmingly begin and end with recognizable dinucleotide signals, and known annotations can supply expected junctions in advance.

This is why the next chapter’s recommendation for RNA-seq differs from its recommendation for DNA. Not because one tool is faster, but because a contiguous aligner is answering the wrong question.

6.8 What to take forward

Alignment converts strings into coordinates, and everything downstream trusts those coordinates without re-deriving them.
Exact alignment does not scale, so every tool seeds and extends. The assumptions in the seeding step determine what each tool is good for.
Error rate, not read length as such, drives the choice between exact-match and sketch-based indexing.
Mapping quality answers a different question from alignment score, and filtering on it is a methodological choice.
A BAM’s properties — sorted, indexed, which reference — are part of its type.

Section 7.1 takes each of these and names what to run.