The matrix a corpus wrote

Every cost model on this ladder has been a rule written down and then measured against. A cost that is not one priced substitutions by the log-odds of a stated random drift; the zero that moves the answer changed what a gap costs at the ends; the edit that reaches back two rows added a transposition at cost one. In each case the numbers came from a rule the essay stated, which is the right discipline for a collection that refuses to quote tables it cannot check.

That essay ended by pointing at the question the discipline had avoided. In every serious application of alignment, the substitution costs are not written down. They are estimated — from a corpus of sequences whose correspondences are already known, by counting how often each pair of symbols is aligned and comparing that count with what chance would produce. The resulting matrix is therefore a statement about the corpus, not about the symbols. This page does that estimation, from corpora whose properties are stated, and measures how much the answers depend on which corpus was used.

How a matrix is fitted #

The procedure is the one used to produce real substitution matrices, reduced to its arithmetic.

Start with a corpus of pairs of sequences, each pair already aligned column by column. Count, over every column of every pair, how often each pair of symbols appears together — $a$ opposite $g$, $c$ opposite $c$, and so on — and how often each symbol appears at all. The count for a pair, divided by the total, estimates the probability $P(x,y)$ that a column holds that pair. The product of the two symbols’ individual frequencies, $p(x)\,p(y)$, estimates the probability that the pair would appear if the two sequences had nothing to do with each other.

The log of the ratio of those two probabilities says how much more or less often than chance the pair is aligned. A pair aligned more often than chance is evidence of correspondence and should be cheap to align; a pair aligned less often than chance should be expensive. So the cost of aligning $x$ against a different symbol $y$ is the negative of that log-odds, in bits, rounded to a whole number and floored at one — a pair aligned far less often than chance is dear, and a pair aligned about as often as chance costs the minimum. Aligning a symbol with itself costs zero, the convention every stated model on this site uses, so a fitted matrix and a written one can be compared entry for entry. A gap character costs three throughout. A pseudocount of one is added to every pair’s count so that a pair never seen in the corpus gets a large finite cost rather than an infinite one.

The corpora are generated so that their properties are known exactly. Each pair is a random sequence and a copy of it that has drifted: at each position the symbol stays with probability stay, and otherwise moves to another symbol, eight times more readily within $\{a, g\}$ or $\{c, t\}$ than across them, with occasional insertions and deletions. A corpus with a stay probability of 0.9 is a collection of close relatives; one of 0.5 is a collection of distant ones.

Two divergences, two matrices #

The two matrices are fitted by the same procedure to the same number of columns and they disagree about every off-diagonal entry.

The near corpus makes substitutions expensive, because in a corpus of close relatives a substitution is rare: every pair of different symbols is aligned far less often than two independent symbols would be paired by chance, so its log-odds is strongly negative and its cost large. Within that, it reproduces the drift’s structure — $a$ against $g$ and $c$ against $t$ cost 2, because the drift moves within those pairs eight times more readily, and the pairs across them cost 4 or 5.

The far corpus makes every substitution cheap, because in a corpus of distant relatives substitutions are common and each pair of different symbols is aligned nearly as often as chance would pair them. The structure is still there in the counts, but after rounding to whole bits there is only room for costs of 1 and 2.

Neither matrix is wrong. Each is the correct estimate of the same model for the corpus it was fitted to. The difference is exactly a model is not a corpus: the log-odds model is one thing and its parameters are a fact about the data the model was shown.

What a bit of cost does to an alignment #

The two matrices matter to an aligner through one comparison above all: a substitution’s cost against a gap’s. A gap character costs three under both. Under the near matrix a transversion costs five, more than a gap character, and a transition costs two, less than one. Under the far matrix every substitution costs one or two, always less than a gap character.

That changes what a minimising alignment prefers in a specific, predictable way. Consider two columns where the strings disagree by a transversion. Under the near matrix, a substitution there costs five, while an insertion and a deletion that realign the strings around it cost six — so the substitution still wins, narrowly, and a small change elsewhere in the pair can tip it. Under the far matrix, the substitution costs two and the two gaps cost six, and nothing tips it. So the near matrix puts many alignments close to a tie between a substitution and a pair of gaps, and the far matrix puts almost none there. An aligner using the near matrix produces more gaps on pairs that are not near relatives, because it has been told that substitutions are rare — which is true of the corpus and false of the test pairs.

This is the same mechanism the zero that moves the answer found for gap costs at the ends of a string: the alignment is decided by comparisons between costs, not by any one cost, and a change that alters one side of such a comparison moves every alignment that was close to it. The parameter plane of the parameter plane has few answers draws those comparisons as boundaries.

The three matrices show the dependence is graded rather than binary. As the corpus gets more divergent, every substitution gets cheaper, and the difference between the cheap and dear substitutions shrinks.

That is the pattern real matrices are organised around. The widely used families of protein substitution matrices are indexed by divergence — one family by the evolutionary distance of the sequences it was estimated from, another by the identity threshold used to cluster them — and choosing a matrix from such a family is choosing an assumption about how far apart the sequences being aligned are. Aligning close relatives with a matrix fitted to distant ones undercharges every substitution, and the reverse overcharges them.

A small corpus writes its own noise into the matrix #

The second source of variation is size. Eight pairs at the same divergence contain 471 aligned columns, and many of the off-diagonal pairs appear only a handful of times in them. A pair seen twice instead of the expected four has its estimated probability halved and gains a bit of cost; a pair seen six times loses one. The small corpus’s matrix has a cheaper entry and a dearer entry than any in the large corpus’s, and both are sampling accidents.

A pseudocount softens this and does not remove it. With a pseudocount of one, a pair never seen still gets a large finite cost — one way a small corpus reaches a cost of 6 — and a pair seen once is treated as twice as frequent as it was. The choice of pseudocount is itself a parameter that moves the matrix, and it matters most exactly when the corpus is small enough to need it.

This is the familiar shape of fitting, and it has the familiar name: a matrix estimated from too little data overfits it.

The arithmetic says how much data is too little. A cost is a rounded logarithm of a ratio of probabilities, so an entry changes by one bit when its estimated probability changes by a factor of about two. A count of $c$ observations has a relative sampling error of about $1/\sqrt{c}$, so an entry is stable to within a bit only once its count is large enough that $1/\sqrt{c}$ is well under a half — a few dozen observations at the least. In eight pairs at a stay of 0.9, about 47 of the 471 columns are substitutions, spread over twelve ordered pairs of symbols and weighted heavily towards the two transitions; a transversion pair is seen perhaps once or twice. Its cost is then a rounded logarithm of almost nothing, and a coin flip in the corpus is a bit in the matrix. It describes the eight pairs well and the drift that produced them poorly. Fitting a class to measurements makes the same argument about fitting growth rates to a few points on this site, where a curve that passes through every measurement is a worse description of the algorithm than one that does not.

Whether the alignments move #

Matrices that differ in their numbers matter only if they differ in their answers. An alignment is a minimum over paths through a table, and many different cost matrices produce the same minimising path — so the question to ask is not how different the matrices look, but how many optimal alignments change when one replaces the other.

The test set is 200 new pairs of forty symbols, drawn at a stay probability of 0.7 — between the near and far corpora — with slightly more insertions and deletions. Each pair is aligned optimally under the matrix from the large near corpus, and then under each other matrix, and an alignment counts as moved if its two aligned rows differ.

More than half of the alignments change between the near and far matrices. Those are not the same alignments with different scores; they are different claims about which symbols correspond — where a gap goes, whether two symbols are a substitution or an insertion followed by a deletion. Anything downstream that reads the correspondences — a variant call, an inferred edit history, a consensus — reads different data depending on the matrix.

The small corpus is the more uncomfortable result. At the same divergence, a matrix fitted to eight pairs moves 79 of the 200 alignments, two fifths of them. The eight pairs were drawn from exactly the drift the large corpus was drawn from, so nothing about the underlying relationship changed. The matrix changed because the estimate was noisy, and the alignments followed the noise.

Forty pairs is ten times the small corpus and a tenth of the large one, and its matrix differs from the large corpus’s in only two entries. Those two entries move 60 alignments — three in ten.

That ratio — two cost entries, thirty per cent of the answers — is the finding the rest of the page rests on. It means the optimal alignments of ordinary test pairs sit close to ties, where a one-bit change in a single substitution cost is enough to make a different path cheapest. An alignment is a discrete answer computed from continuous evidence, and when the evidence is balanced, a small change in the parameters flips the answer completely rather than nudging it.

What is being chosen, when a matrix is chosen #

Three things follow for anyone reading an alignment produced with a fitted matrix.

The matrix is an assumption about divergence. It was fitted to sequences at some distance from each other, and it is the correct model only for sequences at about that distance. A report of an alignment that does not say which matrix, and therefore which divergence, was assumed has left out a parameter that changes more than half the answers.

The matrix carries its corpus’s sampling error. A matrix fitted to a small or unrepresentative corpus has entries that are off by a bit here and there, and a bit is enough to move a substantial fraction of alignments. A corpus that was not generated measured, for this collection’s text structures, how far a real corpus can sit from the generated ones an analysis assumes; a real alignment corpus has the same problem with the additional twist that its pairs were aligned by some earlier method, using some earlier matrix.

An optimal alignment is optimal for a matrix. The phrase “the optimal alignment” names an argmin, and an argmin moves with its objective. A cost that is not one found two models giving two alignments of one pair; this page finds that two estimates of one model do the same thing to a third of a test set. A distance divided by a length found the same fragility in the objective’s form; this is the fragility of its parameters.

What the measurement leaves out #

The corpora are perfectly aligned. Every pair was generated with its true alignment known. A real corpus’s alignments were themselves computed with some matrix, which feeds the matrix’s assumptions back into the counts used to fit the next one. Iterating that loop is how real matrices were refined, and whether it converges to the generating model or to a fixed point of the method is a separate question.

Only substitutions are fitted. The gap cost is held at a fixed value throughout, and gap costs are usually chosen separately, often by trial on a benchmark. How the fitted substitutions and a chosen gap cost interact is exactly the question the next rung takes up.

Four symbols. A protein alphabet has twenty, so a matrix has 190 distinct substitution costs rather than six, each estimated from far fewer columns per entry at the same corpus size. Everything about sampling noise on this page is worse there by a large factor.

Generated corpora. Every corpus here came from one stated drift process, so the only differences between corpora of the same divergence are sampling. A real corpus mixes pairs at many divergences and from many processes, and a matrix fitted to it is an average over them. What the generated collection was right about is this collection’s audit of when generated data predicts real data well, and the answer there was that it predicts the mechanisms and misestimates the magnitudes — which is the right expectation for the numbers on this page too.

The fitted matrix is also, formally, a model of the corpus in the same sense a compressor’s statistics are. The model is the compressor makes that identity explicit for text: a table of symbol probabilities fitted to data is both a description of the data and a code for it, and its cost entries are code lengths. A substitution cost of five bits is a claim that seeing that pair is a five-bit surprise in this corpus — and a surprise is always relative to what was seen before.

Where this ladder goes next: the parameter plane has few answers #

The last plate implies something about the space of cost models that is worth measuring directly. If two entries changing moves thirty per cent of alignments, then as the parameters move continuously the optimal alignment of a single pair must jump between a small number of discrete answers, holding each over some region of the parameters.

That can be drawn. Fix a pair of sequences and sweep two costs — the cost of opening a gap against the cost of extending one, or a substitution against a gap — over a grid of hundreds of settings, colouring each setting by which alignment is optimal there. The prediction from this page is that the plane divides into a handful of regions, that most of the plane is one alignment, and that the boundaries are exactly where the ties live. Under a linear cost model the prediction is sharper still: doubling every cost changes no alignment, so the regions must be wedges through the origin, and one of the two parameters is not really a parameter at all.