Why alignment exists and what each edit represents

You have two sequences in front of you. Maybe two copies of the same gene from two species, maybe a read off a sequencer and the reference genome it came from. They are almost the same and not quite. The question underneath a huge slice of bioinformatics is this: which position over here lines up with which position over there, and what happened in between. That lining-up is called alignment, and this lesson is about what it is really claiming, before we spend the next two lessons on how to compute it fast.

Sequences drift apart one edit at a time

Recall the edit menu from S7.1. DNA changes in exactly three small ways. A substitution swaps one base for another. An insertion adds one or more bases, a deletion removes them, and the two together are indels. Every difference between two related sequences, however tangled it looks, is some pile of those three moves applied over many generations.

Now flip the problem around. In S7.1 you held the original and applied an edit. In the real world you get the two end products and never watched the edits happen. You hold two present-day strings that both descend, you suspect, from one ancestral sequence you will never see. Alignment is the attempt to reconstruct the correspondence between them: to say, position by position, this base here goes with that base there, and right here a base was gained or lost. It is reverse-engineering the edit history from the output alone.

An alignment is a stack of hypotheses

Write the two sequences one above the other and slide them until they register. Where they will not fit cleanly, insert a gap, drawn as a dash, to push the rest back into line. What you get is a stack of columns, and here is the move that makes alignment more than formatting: each column is a claim about a biological event.

A match column has the same residue top and bottom. The most economical story is that this position was conserved, carried down unchanged from the ancestor to both descendants. A mismatch column has two different residues at the same position. The economical story is a substitution: somewhere on one of the two lineages, that base was swapped. A gap column has a residue on one side and a dash on the other. That is an indel: one lineage inserted a base or the other deleted one.

So the three things you can see in an alignment, match, mismatch, and gap, map one to one onto the three things that can happen to a sequence: conservation, substitution, and indel. The visible marks are a readout of invisible history. That is the whole idea. Everything else is machinery for choosing the best such readout.

Notice an honest wrinkle already. A gap in the top row could mean the bottom lineage inserted a base, or the top lineage deleted one. With only two sequences in hand you cannot tell which. The column tells you an indel happened. It does not tell you on which branch, or in which direction. Hold that thought.

Why you cannot just compare position by position

A programmer's first instinct is to skip the gaps: walk both strings together, index by index, and count how many positions differ. That works right up until the first indel, and then it fails completely, for exactly the reason you traced in S7.1. Delete one base near the front of a sequence and every position after it shifts by one. This is the frameshift picture again, only now it is your comparison that gets knocked out of register, not a ribosome's reading frame. From the indel onward, a naive walk compares base 20 against base 21, base 21 against base 22, and reports a wall of mismatches where the two sequences are in fact still identical, just offset.

Alignment exists to undo that offset. By inserting a single gap at the right place, it slides everything after the indel back into correspondence, and the wall of false mismatches collapses back into matches. That is the core job: gaps are not cosmetic, they are the mechanism that re-registers two sequences an indel has knocked out of phase. And re-registering two texts by inserting the smallest sensible set of gaps is a problem you have already trusted a tool with.

Here is the readout made literal. Once the gaps are placed, interpreting the result is trivial. Walking two aligned strings column by column and naming each event is a handful of lines. The aligner's hard job is choosing where the gaps go, which is the next lesson. Read this and predict its output before you trace it.

def read_alignment(top, bottom):
    # top and bottom are two aligned strings of equal length.
    # A dash "-" marks a gap in that row.
    for a, b in zip(top, bottom):
        if a == "-" or b == "-":
            event = "indel: one lineage gained or the other lost a base"
        elif a == b:
            event = "conserved: unchanged since the common ancestor"
        else:
            event = "substitution: a base swapped on one lineage"
        print(a, b, event)


# An aligned pair. The single gap re-registers everything after it.
top    = "ACGTACGT"
bottom = "ACG-ATGT"
read_alignment(top, bottom)

Homology is not similarity

Now the distinction that trips up nearly everyone and that the rest of your bioinformatics career leans on. Similarity is a measurement you can take right now: what fraction of the aligned columns match, or more generally what the alignment scores. It is a property of the two strings in front of you. Homology is a claim about the past: that the two sequences share a common ancestor. One is arithmetic. The other is history.

They are related, and that is exactly why they get confused. High similarity over a long stretch is strong evidence of homology, because the chance of two long sequences matching closely by accident falls off fast as they get longer. But evidence is not identity. Two short sequences can look similar by pure chance, with a four-letter alphabet and not many positions to fill. And the reverse happens too: two sequences can be genuinely homologous yet look barely similar, because so much time and so many substitutions have piled up since they split that the family resemblance is nearly washed out. Similar and homologous are not the same word, and neither one strictly implies the other.

Make it concrete. The widget below runs a real alignment. You will drive its algorithm and scoring choices in the next lesson, so here just read its output. It opens on the classic pair GCATGCU against GATTACA.

Do three things. First, look at the reconstructed alignment rows and, going column by column, name each one as a match, a mismatch, or a gap, and say which event it stands for: conservation, substitution, or indel. Second, edit the two sequences to be nearly identical and watch the columns fill with matches while the score climbs. Third, replace them with two stretches of unrelated-looking bases and notice the crucial thing: the aligner still returns an alignment and a score. It never refuses. That refusal to refuse is precisely why the similarity-is-not-homology warning matters.

matrix[i][j] = max(
  matrix[i-1][j-1] + s(a[i-1], b[j-1]),   // diagonal: align the two residues
  matrix[i-1][j]   + gap,                 // up: a gap in sequence b
  matrix[i][j-1]   + gap                  // left: a gap in sequence a
)
s(x, y) = match when x == y, else mismatch

// A thin sketch of needlemanWunsch (global). Local alignment is the same
// loop with a 0 floor and a traceback that starts at the largest cell.
function fill(a: string, b: string, match: number, mismatch: number, gap: number) {
  const m = a.length;
  const n = b.length;
  const matrix: number[][] = [];
  for (let i = 0; i <= m; i++) matrix.push(new Array(n + 1).fill(0));

  // Seed the borders: a run of k gaps costs k * gap. (Local leaves these at 0.)
  for (let i = 1; i <= m; i++) matrix[i][0] = i * gap;
  for (let j = 1; j <= n; j++) matrix[0][j] = j * gap;

  // Fill. Row i indexes a, column j indexes b.
  for (let i = 1; i <= m; i++) {
    for (let j = 1; j <= n; j++) {
      const s = a[i - 1] === b[j - 1] ? match : mismatch;
      const diag = matrix[i - 1][j - 1] + s;   // align a[i-1] with b[j-1]
      const up = matrix[i - 1][j] + gap;        // a gap in b
      const left = matrix[i][j - 1] + gap;      // a gap in a
      matrix[i][j] = Math.max(diag, up, left);
      // Smith-Waterman (local) instead: Math.max(0, diag, up, left)
    }
  }

  // The global score is the bottom-right corner.
  return matrix[m][n];
}