What two letters can carry

The copy rate \(\kappa\), by doing

The coin of Part 1, section 1 fixes the match rate inside a shared domain. Heads, with probability \(\kappa\), copies the query's class into the target. Tails draws the target's class from its own composition, and a tails position still matches by chance, with probability \(\Pr(\text{match} \mid \text{unrelated})\). Together:

(1a)\[\begin{aligned}\Pr(\text{match} \mid \text{shared domain}) &= \kappa + \bigl((1 - \kappa) \cdot \Pr(\text{match} \mid \text{unrelated})\bigr) \\ &= \Pr(\text{match} \mid \text{unrelated}) \\ &\quad + \Bigl(\kappa \cdot \bigl(1 - \Pr(\text{match} \mid \text{unrelated})\bigr)\Bigr)\end{aligned}\]

Solving Equation (1a) for \(\kappa\) gives Equation (1b):

(1b)\[\kappa = \frac{\Pr(\text{match} \mid \text{shared domain}) - \Pr(\text{match} \mid \text{unrelated})}{1 - \Pr(\text{match} \mid \text{unrelated})}\]

Equation (1b) is Cohen's \(\kappa\), and it is the copy rate: the match rate above chance, as a fraction of the room above chance. Measured on a pair built with the coin, it gives the coin's \(\kappa\) back; Figure 1, just below, lets you build one and check. 0.46 is the value measured for hp_thomas_dill2 in aligned Pfam pairs at 20–30% identity (analysis notebook 230). The expected score per step inside a shared domain is \((1 + C)\,\Pr(\text{match} \mid \text{shared domain}) - C\); the note under Figure 1 of Part 1 prints it for the walk.

a match: same class in both sequences a mismatch outlined: the coin copied this position (known only in the simulation) hatched bar: matches that come for free, the chance match rate bar: matches above chance

Figure 1. Top: the real pair, residues, their classes in the picked alphabet, and a match or mismatch at each position. Middle: its match rate as a bar from 0 to 1; the hatched part up to the chance rate comes free, and κ is the fraction of the room above chance that the pair fills. Bottom: a simulated pair on the same chance rate, with the coin's copied positions outlined, and κ recovered from its counts. The simulated pair is a picture of the definition and feeds nothing else. λ_pair (Part 1, section 3) is solved from the class shares alone and describes unrelated sequences, κ = 0. K and r_database (Part 1, section 5) are fitted on real database sequences, not on this.

κ, the copy rate of the simulated pair (Cohen's κ) 0.46

alphabet the pair is recoded in; the chance match rate follows from its class shares 0.50

First a real pair: human BCL-2 residues 130–166 against C. elegans CED-9 residues 154–190, the BH1 region, 5 of 37 residues identical. It is recoded in the alphabet you pick. The hp tables differ on C, G, P, W and Y, and this window has all five, so the count moves with the table. Each position either matches or does not, and that is all a real pair can tell you: a match that was conserved and a match that happened by chance look the same. What you can do is count. The match rate is matches over positions. The chance match rate comes from the two sequences' own class shares, the way Part 1, section 2 explains. κ is the gap between them as a fraction of the room above chance, Equation (1b).

Then a simulated pair on the same chance rate, built with the coin at the copy rate you set. Here, and only here, the copied positions are known, because the simulation flipped the coin; they are outlined. κ recovered from the simulated counts lands back on the copy rate, up to the luck of the draw.

Equation (1b) for the real pair and for the simulated pair:

\[\kappa\]

real pair: matches of 37 positions—

real pair: chance match rate from its class shares—

real pair: κ, Eq. (1b)—

simulated: copied / lucky match / mismatch—

simulated: match rate—

simulated: κ recovered, against the copy rate set—

The same window as evidence, against the database and query length set beside Figure 2 (Equations (2a) and (2b) below, at this pair's own κ and class shares):

bits per position this window carries at its own κ (mutual information, Eq. 2a)—

bits in all 37 positions—

bits that E = 1 costs, log₂(K · m · n)—

the window scored +1 / −C: S, and E for that S—

—

Measured, not modelled: searched against the 19,732 canonical human proteins with every kmerseek alphabet at k set so one seed carries 16 to 44 bits, and C at C_best, the best E-value this pair reached was 11,200 in two letters (hp_lehninger2, k = 19) and 1,950 in any alphabet (wwmj5, k = 9); BCL-2 ranked 1,352nd of 5,153 human proteins hit under the E-value, and 213th of 18,064 under the mean IDF of the region's k-mers (polarity4, k = 9). Analysis notebook 241, 2024-kmerseek-analysis PR #46. The E-values are far above 1, as the readout above predicts. The ranks say what a window this short can still do: put the partner in the top 1 to 25% of the proteome, not at the top.

Bits per position: how many bits two letters carry, and how long a twilight-zone match must be

two letters at this κ every other alphabet at its measured κ, or, where none is measured, at a κ estimated from its number of classes on the line between the measured alphabets (names on the right; bold = measured) BLOSUM45, 0.38 bits per position (the 20-letter matrix for 20–30% identity) bits that E = 1 costs, log₂(K m n) Pfam domain length, 100–200 residues

Hover a name on the right to read its bits per position; click one to make it the alphabet of the whole page.

Figure 2. Equation (2a) and (2b) drawn. A 150-residue domain at 20–30% identity carries about 24 bits in two letters, against the 24 to 35 bits that E = 1 costs on the three databases (log₂(K m n) at m = 240, the starting value of the slider); the same domain carries about 57 bits with BLOSUM45. The two lines cross the E = 1 line at the lengths marked. The thin lines are every other alphabet of Part 1, section 2, each at its own measured κ where one exists and otherwise at a κ estimated from its number of classes (the line between the four measured alphabets, against log classes), so the fan shows the whole trade at once: more classes, more bits per position, and none reaches the full matrix. Drag the coral line to set κ: its slope is the bits per position, and κ is found from that. The estimate uses one κ for every position; the text says in which sense it is a floor and in which a ceiling.

\(\kappa\) is the match rate above chance between aligned sequences that share a domain, 0.46 for two letters at 20–30% identity. It is not the same thing as information. The information an aligned position carries is the mutual information between the query class and the target class, and that is what any scoring system, at best, collects. The page models a pair that shares a domain with the coin of Part 1, section 1. Heads (probability \(\kappa\)) copies the query's class into the target, a conserved position. Tails draws the target's class from its own composition regardless of the query, a drifted one. \(\Pr(\text{match} \mid \text{shared domain})\) then obeys Equation (1a), and the bits follow from the joint distribution of the two classes, Equation (2a):

(2a)\[\text{bits per position} = \sum_{i,j \,\in\, \text{classes}} \Pr(i,j)\,\log_2 \frac{\Pr(i,j)}{\Pr(i)\,\Pr(j)}\]

Divide the bits the search space costs, \(\log_2(K \cdot m \cdot n / E)\), by the bits per position, and you get how long an ungapped stretch has to be, Equation (2b):

(2b)\[\text{positions needed} = \frac{\log_2\bigl(K\,m\,n / E\bigr)}{\text{bits per position}}\]

What the numbers say: at \(\kappa = 0.46\), the value measured for hp_thomas_dill2, two letters carry about 0.16 bits per position, against 0.38 for BLOSUM45 in 20 letters. \(E = 1\) costs \(\log_2(K\,m\,n)\) bits: about 24 on SCOPe40, 25 on one proteome, 30 on Swiss-Prot and 35 on UniRef50 for a 240-residue query. So an ungapped stretch needs about 150 positions on SCOPe40, 160 on one proteome, 190 on Swiss-Prot and 220 on UniRef50; the readouts under Figure 2 recompute these for the sliders. A 200-residue domain at ordinary conservation gets through on SCOPe40 (\(E \approx 0.004\)) and just about on Swiss-Prot (\(E \approx 0.3\)). A 150-residue one sits at the threshold on SCOPe40 (\(E \approx 1\)) and is lost on the proteomes (\(E \approx 60\) on Swiss-Prot, 3,000 on UniRef50). Pfam domains run 100 to 200 residues, and gaps break the stretch.

So on a database the size of Swiss-Prot or UniRef50 most of the twilight zone is out of reach for a two-letter search. On SCOPe40, or on one proteome, which costs about the same, the longer half of it is in reach, but only for a region that stays ungapped for 150 residues or more. A pair whose conserved part is one short block is not: the BH1 window of Figure 1 is 37 residues and carries — bits, and no proteome-sized database costs less than 24. Out of reach means the E-value cannot single the pair out from the chance regions that share its seed. The seed itself is found: the —-residue class-identical run between BCL-2 and CED-9 under Figure 1 is in the index, together with about — unrelated places in — that carry the same class string (n / 2^—). Against a smaller search space, one family instead of a proteome, the same evidence can be enough, because \(E\) scales with \(K \, m \, n\).

Ranking the hits by how rare their shared k-mers are (mean IDF, tf-idf, enrichment) does not change this. IDF is \(\log(n / \text{count of the k-mer})\), the same search-space term \(E\) already contains. In two letters every 19-mer is about equally common, so a "rare" two-letter region is one of unusual composition, not one of shared ancestry. Re-weighting the same shared positions cannot add bits the pair does not carry. Measured on 687 Pfam pairs, the E-value ranks true regions best (average precision 0.82) and mean IDF worst (0.22). On the BCL-2/CED-9 pair against the human proteome, mean IDF moved BCL-2 from 1,352nd to 213th, still not first (2024-kmerseek-analysis, PR #44 and notebook 241).

The estimate is a ceiling in one sense and a floor in another.

Too high, with respect to scoring: \(+1/-C\) collects only — bits per position at these settings, under the — bits the mutual information allows, because a flat penalty cannot weight a rare-class match above a common one.
Too low, with respect to conservation varying along the domain. The coin uses one \(\kappa\) for every position. Real domains have a conserved core and drifted loops, and mutual information is convex in \(\kappa\), so a mix of strongly and weakly conserved positions carries more than the single-\(\kappa\) figure. Pairs with a conserved hydrophobic core or a strong amphipathic pattern carry more still, and those are the pairs a search finds.
Too low, with respect to what is scored: one ungapped region. Summing several regions in consistent order (section 5), a matrix over classes in place of \(+1/-C\), or a profile over the class string all collect more of the available bits.

Separately from bits, the 12-mer seed is its own filter: the median longest class-identical run at 20–30% identity is 12 residues, so about half of true pairs never produce a seed however many bits they carry. Switch the alphabet in the bar at the top to see the bits per position grow with the number of classes.

Can two letters see what twenty cannot?

the coin model at each alphabet's measured κ: bits per aligned position on the left, expected longest class-identical run in a 150-residue domain on the right BLOSUM45, a full 20-letter matrix k = 12, the seed

Figure 3. The trade the alphabet makes. Left: bits per aligned position rise with the number of classes, and no reduced alphabet reaches the full matrix. Right: the longest run of class-identical positions falls with the number of classes, and only two letters reach the 12-mer seed at 20–30% identity. Both panels use the Swiss-Prot average composition and the κ measured for each alphabet in aligned Pfam pairs at 20–30% identity.

No: per aligned position, two letters cannot see what twenty cannot. This is a theorem, not a measurement. The class of a residue is a function of the residue, so merging letters can only lose mutual information, never create it. The twilight zone is named for 20-letter identity: below about 30% identity, alignment in the full alphabet stops being able to tell relatives from chance. The hope is that a two-letter string, where the same pairs still match at 70% of positions, sees through it; the theorem rules that out. Figure 3 puts numbers on it with the coin model at each alphabet's own measured \(\kappa\). Two letters carry 0.16 bits per position, four 0.25, and twenty 0.28 under the coin model, which treats every substitution as a miss. BLOSUM45, which scores twenty letters with a full matrix, collects 0.38. A 70% class match rate against a 50% chance rate is worth less than a 25% residue match rate against a 6% chance rate.

What a small alphabet buys is not bits but seeds. kmerseek finds a region only through an exact shared 12-mer, and the expected longest class-identical run in a 150-residue domain at 20–30% identity is about 12 in two letters, 9 in four, 6 in five and 3 in twenty. Exact 12-mers exist at that divergence only in the two-letter string. So a two-letter k-mer index can find pairs that a twenty-letter k-mer index cannot find at all, and it does so without alignment. But each pair it finds carries fewer bits than a matrix-scored alignment of the same pair would, and the E-value says so. That sets the fair claim for kmerseek. Against another k-mer index (a twenty-letter k-mer search, which needs exact twenty-letter k-mers), two letters reach lower identity, at k-mer speed. Against BLAST, which aligns with a full substitution matrix and needs no exact k-mer, two letters do not see more: they see less per position, and only faster.

Short domains and motifs, 20–40 residues

the alphabet chosen in Part 1, section 2, scored +1 / −C at the class match rate you set twenty letters under the coin model at 20–30% identity, every substitution a miss BLOSUM45, a full matrix what the search space costs, ln(K · m · n)

Figure 4. A motif of the length you set, at the class match rate you set. Left: the evidence it carries in each alphabet against the cost of the search space; a bar that does not reach the dashed line has E above 1. Right: whether a class-identical run of k residues, the seed kmerseek needs (k = 12; the slider tries others), is likely inside it, and how many such k-mers the database holds by chance. The expected longest run and the probability of a run of at least k are computed without approximation, not from the large-length formula, which goes wrong when the motif has fewer than one expected mismatch. Press Grow to watch the motif lengthen.

motif length, residues 30

class match rate inside the motif (1 = every position) 0.75

k, seed length for the seeding question (kmerseek's calibration and every E-value here use k = 12) 12

Bits per position set a floor on length, and a 20–40 residue motif is under it in every alphabet. At twilight-zone conservation (class match rate 0.75), 30 positions carry about — bits in two letters, 8 in twenty letters under the coin model and 11 with BLOSUM45, against the — bits the search space costs on —. Set the class match rate to 1 to see what a perfect class match is worth: 30 residues give E ≈ —, 40 residues E ≈ —.

Pairwise search of any kind, in any alphabet, is the wrong tool for a short motif at that identity, which is why Pfam finds such domains with profiles built from many sequences. What a two-letter index can still do is find a short conserved motif exactly: a 12-residue class-identical run is a seed. But a two-letter 12-mer occurs by chance about n / 2¹² times in a database of n residues, 600 times in SCOPe40 (the figure recomputes this for the alphabet chosen in Part 1, section 2), so the seed alone is findable, not assessable.

Equations (1a), (2a) and (2b) at the κ, shares and database set beside Figure 2:

\[\text{bits per position}\]

shared-domain match rate at this κ—

bits per aligned position, this alphabet (mutual information)—

bits per position that +1 / −C collects—

bits needed for E = 1 (E = 0.001)—

aligned positions needed, no indel—

an exact class-identical run alone—

Several regions in a row: what summing adds, and what each extra region costs

BCL-2 and CED-9 share four short conserved blocks, BH4, BH3, BH1 and BH2, in the same order along both proteins, and kmerseek reports each as its own region. Adding their scores is the idea behind the sum statistics of Karlin and Altschul (1993). A second region only counts as evidence when it sits where a shared domain would put it: after the first, in both sequences, within some gap allowance \(g\). Saying where it sits costs bits, about \(\log_2 g\) for its place in the query and \(\log_2 g\) for its place in the target. So \(r\) regions of \(b\) bits each are worth about Equation (3). This is the page's own rough count; the exact statistic is in the reference.

(3)\[\text{bits from } r \text{ regions} \approx (r \cdot b) - \bigl((r - 1) \cdot 2 \cdot \log_2 g\bigr)\]

\(r\): regions in consistent order. \(b\): bits one region carries (its score times λ_pair times r_database, in bits). \(g\): the largest gap allowed between consecutive regions, in residues, in either sequence.

outline: r × b, the bits the regions carry before placement filled: what is left after paying (r − 1) × 2 log₂ g to place them bits that E = 1 costs, log₂(K m n), on the database set beside Figure 2

—

Figure 5. Equation (3) drawn for one to eight regions. Each region carries the bits set on the slider; the starting value is what the BH1 window of Figure 1 scores at +1 / −C. The filled bar rises with r only while one region carries more than the 2 log₂ g bits its placement costs; below that, adding regions loses evidence. The dashed line is the same E = 1 cost as in Figure 2.

b, bits per region 9

g, gap allowed between regions (residues) 50

Four regions at the BH1 window's — bits, with a gap allowance of — residues, sum to — bits after placement, against the — bits E = 1 costs on —. Summing helps only when each region is long enough to pay for itself, about — bits at this gap. — Two levers add evidence without adding positions. One is on the query side. Search several sequences that carry the same domain, the query's orthologs or paralogs, the way a Pfam profile is built from many sequences, and keep the targets every one of them finds at the same place. Each query is a separate draw of the same conserved core, and a target hit by all of them shares the domain, whatever its ancestry. The other is the alphabet. Every alphabet merges different residues, so the chance hits at the top of one alphabet's list differ from those at the top of another, while a true match stays high in many. A rank combined across alphabets can then sort the partner above the chance hits even when no single E-value can. Neither raises the bits a pair carries above the 20-letter matrix; both change the rank, not E. A mimic, which shares a surface patch or a short motif and not a domain, is the case of section 4: findable as a class-identical run, not assessable by any pairwise E-value.

The best mismatch penalty \(C_{\text{best}}\) shrinks as the alphabet grows

The scoring system that finds a shared domain most efficiently gives each pair of classes the log of how often they meet in a shared domain over how often they meet by chance, \(\ln\bigl(\Pr(i, j \mid \text{shared domain}) / \Pr(i, j \mid \text{unrelated})\bigr)\); this is the result behind BLOSUM (Altschul 1991). In the coin model of \(\kappa\) above, two different classes meet in a shared domain only on a tails flip. So a mismatch scores \(\ln(1 - \kappa)\) whatever the two classes are. A match of class \(c\) scores \(\ln(1 + \kappa / f(c) - \kappa)\), where \(f(c)\) is the share of class \(c\): rarer classes, stronger matches. kmerseek scores +1 and \(-C\), so the \(C\) that copies this system is the ratio of the two:

(4)\[C_{\text{best}} = \frac{-\ln(1 - \kappa)}{\ln\bigl(1 + \kappa / f(c) - \kappa\bigr)}\]

C_best against the number of equal classes, at the κ of the slider an alphabet with a measured κ, at its own κ (bold in the list) an alphabet with no measured κ, at a κ estimated from its number of classes (the line between the measured alphabets) C_min = p / (1 − p) with equal shares: below it a chance position scores zero or more, no λ exists, and no E-value can be computed× an alphabet whose C_best at equal shares is under its own C_min at the Swiss-Prot shares: kmerseek at that C has no λ for it the C of this page (2 to start, kmerseek's setting); drag it or the slider and every figure on the page follows

Hover a dot or a name to read it off; click one to make it the alphabet of the whole page.

Figure 6. Equation (4) with equal class shares, f(c) = 1 / classes. The curve is one κ for every alphabet size. Every alphabet of Part 1, section 2 is a dot at its number of classes: a filled dot uses the κ measured for that alphabet in aligned Pfam pairs at 20–30% identity, which itself falls as the alphabet grows; an open dot has no measured κ and uses a κ estimated from its number of classes, on the line between the four measured alphabets against log classes. Alphabets with the same number of classes share a dot and are listed one per row. C = 2 sits above the two-letter value of 1.6 and well above every larger alphabet's. The grey line is the floor: with equal shares a chance position scores +1 with probability 1 / classes and −C otherwise, and that has to average below zero, so C must exceed 1 / (classes − 1). Alphabets with unequal shares have a higher floor, p / (1 − p) with p the chance match rate at the Swiss-Prot shares, and a red × marks those whose C_best falls under it. Drag the coral line to change C for the whole page, or drag the purple curve to change the κ it is drawn at.

κ, the copy rate, for the curve 0.46

C, the mismatch penalty for the whole page (kmerseek: 2; the measured search used C_best) 2

With two letters a match is a coin flip that came up right, weak evidence, so a mismatch must count for a lot beside it. With 20 letters a match is rare by chance and strong evidence, while a mismatch is what most positions of a true pair at 20–30% identity look like, so it should cost little. Any \(C\) gives valid E-values once \(K\) and \(r_{\text{database}}\) are recalibrated at that \(C\); they were fitted at \(C = 2\). \(C\) sets which regions score high, not whether the statistics hold, and \(X\) should scale with it. Above two letters a flat \(-C\) also ignores that some class pairs are near-misses. That is where a matrix over classes would beat any single \(C\).

Equation (4) has a floor. The expected score of one chance position is \(p - C \cdot (1 - p)\), with \(p = \Pr(\text{match} \mid \text{unrelated})\), and Karlin–Altschul statistics need it below zero, so \(C > p / (1 - p)\). With equal shares \(p = 1 / \text{classes}\) and the floor is \(1 / (\text{classes} - 1)\), the grey line. Real shares are unequal, \(p\) is the sum of the squared shares and is larger, and the floor rises. In the measured search (notebook 241) that is what happened for hp_lehninger_hpc3, gbmr4 and gbmr7: at the \(C_{\text{best}}\) derived with equal shares a chance position scored zero or more, no \(\lambda\) existed, and no region of theirs got an E-value at any \(k\). The fix is to derive \(C_{\text{best}}\) from the measured shares, so the note under Figure 6 prints it at the shares of the chosen alphabet.

Glossary

Alphabetical: Latin letters first, then Greek by name, then words. Symbols defined in Part 1 are marked so.

bits per position: the mutual information between the query class and the target class at one aligned position, Equation (2a).
\(C\): mismatch penalty, --extend-mismatch-penalty; 2 unless the slider in section 6 is moved. The measured search of notebook 241 used \(C_{\text{best}}\) for each alphabet. \(C_{\text{best}}\): the penalty that copies the log-odds scoring, Equation (4). \(C_{\text{min}}\): the penalty below which a chance position scores zero or more and no \(\lambda\) exists, \(p / (1 - p)\); section 6.
\(E\): expected number of chance regions scoring at least \(S\), Part 1, Equation (6a). Smaller is better.
\(f(c)\): the share of class \(c\) in a sequence; \(h\) and \(p\) for two classes.
\(K\): the per-database constant in front of \(E\), measured when the index is built; Part 1, section 5.
\(k\): seed length; a region is found only if it contains an exact shared k-mer. 12 on this page; 16 to 24 for two letters in the measured search of notebook 241, where k was set so one seed carries 16 to 44 bits in every alphabet.
\(m\): query length in residues. \(n\): database residues.
\(\Pr(\text{match} \mid \text{shared domain})\): fraction of aligned positions, between two sequences that carry the same domain, where the classes match; Equation (1a).
\(\Pr(\text{match} \mid \text{unrelated})\), the chance match rate: how often one query position and one target position, drawn independently, carry the same class; Part 1, Equation (2b).
\(r_{\text{database}}\): the share of a region's score that counts as evidence, one number per database; Part 1, section 5.
\(r\), \(b\), \(g\): in Equation (3), the number of regions summed in consistent order, the bits one region carries, and the largest gap allowed between two of them, in residues.
\(S\): region score, matches minus \(C\) times mismatches; Part 1, Equation (1).
\(\kappa\) (kappa), the copy rate, Cohen's \(\kappa\): the match rate above chance as a fraction of the room above chance, Equation (1b).
\(\lambda_{\text{pair}}\) (lambda): nats per unit of score for this pair; Part 1, Equation (3).
alphabet: the classes the 20 amino acids are merged into before anything is scored (--alphabet); chosen in the bar at the top.
nat: the natural-log unit of evidence; one nat is \(1/\ln 2 = 1.44\) bits.