Part 1 set out how a k-mer match earns a score and an E-value, and Part 2 worked out how much evidence a reduced alphabet can carry. This page measures it: when three proteins are searched against all 19_732 canonical human proteins with each of kmerseek's 19 alphabets, what does each alphabet see? Every dot below is a real search. Pick an alphabet, a test case and a ranking metric; every figure and every number on the page answers to those three choices. Sections 7 and 8 then ask whether the best ranks beat what a random query gets, and whether a model trained on known Pfam domains combines the alphabets better.
Part 1, how a region gets an E-value · Part 2, what two letters can carry · results notebook 241 in 2024-kmerseek-analysis #46, notebooks 242 and 243 in #48
Each query is searched against the same database, the 19_732 canonical human proteins of GENCODE v49 (11_347_201 residues). The two proteins we are looking for, BCL2 and CD47, are both in it, so the search that finds BHF's partners is the search that has to find them.
For every alphabet the seed length k is set so one seed carries 16 to 44 bits, where bits per position is the entropy of the alphabet's class shares measured on this proteome. That makes a 2-letter k=24 seed and a 20-letter k=6 seed the same amount of information, and it is why the x axis everywhere is bits, not k. Each index was built with region extension at the alphabet's own mismatch penalty and a Karlin-Altschul fit; every metric kmerseek writes was kept for every region, with no threshold.
E-value: matches this good expected by chance in the whole search, from the fitted lambda and K. Mean IDF: how rare the region's k-mers are in the database, averaged. tf-idf: those rarities summed over the region. Enrichment: shared k-mers divided by the number expected by chance. Poisson p-value: the chance of at least this many shared k-mers under a Poisson model.
A target is one human protein, scored by its best region under the metric. The partner's rank is 1 plus the number of targets strictly better; the number tied with it is shown separately, so rank 1 with 40 tied means "tied for first with 39 others", not found.
Hover any arm for its numbers. The bits axis is the same on every figure of this page.
One line per metric: where the partner ranks at each arm. The thick pale line under them is last place, the number of human proteins hit at that arm, so a metric can never sit below it. A metric with no mark at an arm either lost the partner there or, for the E-value, has no Karlin-Altschul fit at that arm.
Each cell is the best rank over all k for that alphabet and metric; the best cell in each row is bold. Click a row to select the alphabet.
with m the query length, n the residues in the database and K taken as 1. The fitted K is smaller, which lowers the bar by a few bits; the gap does not close.
Left: for every arm of the selected alphabet, the share of matched regions whose E-value is infinite. Right: identity. The chance identity of two unrelated sequences under this alphabet is p; a region whose own identity is above C / (1 + C) gets lambda 0 and no E-value.
Two things follow. First, the limit removes the best-matching regions first, because it is the regions that match well that cross it. Second, the kappa-optimal penalty was derived assuming equal class shares (1 / classes), but the real chance identity p is the sum of the squared shares, and for gbmr7, hp_lehninger_hpc3 and gbmr4 it is high enough that a chance position scores at or above zero. There no lambda exists at any k. Deriving C from the measured shares fixes that; the identity limit needs a change to how lambda is solved.
A best rank of 213 out of 19_732 sounds good, but it is the best of 152 arms and five metrics. To see what that is worth, 300 random human proteins of about the query's length (within 25%) were searched in the query's place, with every arm, against the same database. Bcl-2 family proteins and CD47 were left out of the random set. If random proteins rank BCL2 as well as Ced9 does, the rank says nothing about Ced9.
For "best single arm", each query keeps the partner's best rank over every arm of the 19 alphabets, the same choice Ced9 and P66 get. For "all 19 alphabets combined", each alphabet contributes the partner's rank at its lowest-information arm. The 19 ranks are combined by their geometric mean: divided by 19_732, multiplied, and the 19th root of the product scaled back to a rank. A query that does not hit the partner at an arm gets the middle of the tied bottom ranks.
Ranked by mean IDF. A combination's score is the geometric mean of the partner's ranks in its alphabets. Every random protein picks its own best combination, so it gets the same head start as the real query.
Instead of picking alphabets by hand, a logistic regression (a model that weighs each feature and turns the sum into a probability) learned how to combine them. It was trained on 998 human proteins with a Pfam domain, each searched against all the others with every one of the 152 arms. A match between two proteins counts as correct when both matched stretches fall in the same Pfam family, each covering at least 20% of the stretch. The model was trained on half of the Pfam families and tested on the other half, which it never saw. The Bcl-2 family is only in the test half, so Ced9 → BCL2 is a test the model was not fitted to.
Matches from different arms that lie on the same diagonal (within 10 residues) and overlap on the query are pooled into one candidate match. For each of the 19 alphabets the model sees five values: whether any of its arms found the match, the highest seed information at which one still did, the best mean IDF divided by ln(number of proteins in the database), the best bit score and the most shared k-mers. It also sees the length of the matched stretch, the lengths of both proteins, and how many alphabets and arms found the match. Last, it counts on each side the membrane-helix-like windows: 19 residues with at least 14 of A, I, L, M, F, V, W. Each arm was searched with the same setting, extended or exact regions, as in the human search above, so a feature means the same thing in training and in use.
Precision: the share of matches kept that are correct. Recall: the share of all correct matches that are kept. Average precision is the area under a curve. Recall is on a log scale because every curve has fallen close to chance by 30% recall.
The model was then applied to the searches of Ced9, P66 and BHF against the 19_732 human proteins in the sections above. A human protein's score is the probability of its best candidate match.
Click a row to select it. Bits per position is the entropy of the class shares on the human proteome; p is the chance that two unrelated positions fall in the same class, the sum of the squared shares.