How a kmerseek region gets an E-value

• \(X\), the give-up margin (--extend-xdrop, default 8). The walk remembers the best score it has reached. As soon as the current score is \(X\) points below that best, the walk stops and the region ends where the best was. With \(C = 2\), four differing positions in a row cost 8 and end it; a mix of matches and mismatches can carry it further. A larger \(X\) lets the walk cross a short patch of mismatches inside a real shared domain, and also lets it wander further into unrelated sequence before it gives up.
• Shared-domain stretch: a knob of this model, not of kmerseek. It is the number of residues on each side of the seed where the target carries the same domain as the query, so the classes match more often than chance; beyond it the two sequences are unrelated and match only by chance. In a real search this is however far the shared domain extends past the seed.
• \(\kappa\) inside that stretch, the copy rate. The model builds the target from the query one position at a time, flipping a weighted coin at each position. Heads, with probability \(\kappa\): the target gets the class the query has there, as a conserved position would. Tails, with probability \(1 - \kappa\): the query is ignored and the target's class is drawn from the target's own composition, as at a position that has drifted, or between unrelated sequences. So \(\kappa\) is the fraction of positions matching above chance in the chosen alphabet; the rest match only as often as chance allows. At \(\kappa = 0\) the stretch is no different from unrelated sequence. 0.46 is the value measured for hp_thomas_dill2 in aligned Pfam pairs at 20–30% identity (analysis notebook 230). Part 2 shows that this copy rate is Cohen's \(\kappa\) and measures it on a real pair.

Fixed in the walk: \(k = 12\) for the seed, \(C = 2\), and the two compositions, which section 2 lets you set. Draw another pair redraws the random sequences with the same settings.

The shares are counted from the two sequences of the pair: not from the database, and not from a uniform assumption. The sliders stand in for those counts. What the theory assumes is only that each position is independent of its neighbours; that is the part real sequences break, and \(r_{\text{database}}\) in section 5 measures by how much.

The tables are kmerseek's own, from src/rust/alphabets.rs, clusters in the order their classes are numbered; the hp alphabets write their classes as h and p (and c for cysteine in hp_lehninger_hpc3), every other alphabet writes a class as its first residue in lowercase. A F I L M V is hydrophobic and D E H K N Q R S T polar in every hp table; they disagree only about C, G, P, W and Y.

If the average were zero or positive, chance alone would keep the score climbing, a high score would be the normal outcome, and it would mean nothing. A pair of unrelated sequences that already match at more than two positions in three by luck cannot be assessed at all, whatever the score. Whether one region is significant is a different question, answered by \(E\) itself in section 4. Section 3 turns the same requirement into the value of \(\lambda_{\text{pair}}\): the expected score is the slope of the curve in Figure 3 at \(\lambda = 0\), and only a negative slope lets the curve come back to zero.

At \(C = 2\) the positive root has a closed form, Equation (5), written with \(u\) as shorthand for \(\Pr(\text{match} \mid \text{unrelated})\) in this equation and its derivation only:

Four steps, all at \(C = 2\) and with \(u = \Pr(\text{match} \mid \text{unrelated})\). First, put \(C = 2\) into Equation (4):

Second, write \(y = e^{\lambda_{\text{pair}}}\), so \(e^{-2\lambda_{\text{pair}}} = 1 / y^{2}\) and \(y > 0\). Multiply both sides by \(y^{2}\) and move everything to the left. The equation in \(\lambda_{\text{pair}}\) becomes a cubic in \(y\):

Third, \(y = 1\) (that is, \(\lambda_{\text{pair}} = 0\)) is always a root: \(u - 1 + (1 - u) = 0\). That is the root the page keeps saying to ignore. Dividing the cubic by \((y - 1)\) removes it and leaves a quadratic; multiply the right side out to check it gives the cubic back:

Nothing is divided by zero here. The line above is a statement that two expressions are equal for every \(y\), and checking it needs multiplication only. A product is zero when one of its factors is zero, so the cubic is zero when \(y = 1\) or when the quadratic is zero, and the root we want comes from the second case, where \(y\) is not 1 and \(y - 1\) is not 0. The quadratic at \(y = 1\) is \(3 \cdot u - 2\), which is 0 only at \(u = 2/3\): there the root we want has met the one we threw away, and \(\lambda_{\text{pair}} = 0\).

Fourth, solve the quadratic \(u \cdot y^{2} - (1 - u) \cdot y - (1 - u) = 0\) with the usual formula. Under the square root, \((1 - u)^{2} + 4 \cdot u \cdot (1 - u)\) has the common factor \((1 - u)\), and \((1 - u) + 4 \cdot u = 1 + 3 \cdot u\). Taking the \(+\) sign gives the larger root, which is the positive one because the other has a negative numerator:

\(\lambda_{\text{pair}} = \ln y\), which is Equation (5). Check against the boundary of section 2: at \(u = 2/3\) the numerator is \(1/3 + \sqrt{(1/3) \cdot 3} = 4/3\) and the denominator is \(4/3\), so \(y = 1\) and \(\lambda_{\text{pair}} = 0\), which is where the expected score per position reaches zero. Below \(2/3\) the fraction is above 1 and \(\lambda_{\text{pair}}\) is positive. A number to hold on to: at \(u = 1/2\), \(y = 1/2 + \sqrt{5/4} = 1.618\) and \(\lambda_{\text{pair}} = 0.481\) nats per score point; putting \(y = 1.618\) back into the first line gives \(0.809 + 0.191 = 1\). For any other \(C\) the same substitution gives a polynomial of degree \(C + 1\) in \(y\) (when \(C\) is a whole number), which is why the code bisects instead.

Solving \(\lambda\) per pair is in the spirit of BLAST's composition-based statistics (Schäffer et al. 2001), which adjust per query. kmerseek's code and its CSV columns write \(a\) for what this page writes \(\Pr(\text{match} \mid \text{unrelated})\).

\(r_{\text{database}}\): the database factor, the share of a region's score that counts as evidence. The E-value judges a region with score \(S\) as if it had scored \(r_{\text{database}} \cdot S\) between sequences whose positions were independent. On Swiss-Prot, \(r_{\text{database}} = 0.865\): a region scoring 40 is judged like one scoring 35. The 5 points taken off are matches that real proteins hand out for free: a helix, a strand or a membrane stretch keeps one class for several residues in a row, in both sequences, so neighbouring positions match together without the two proteins being related. A plain number with no unit, one per database, measured on that database's own sequences when its index is built (section 5): 1.00 would mean positions are independent; SCOPe40 domains give 0.947, full-length proteins 0.83–0.87. In the code it is ka_r_database. Figure 10 draws the discount and Figure 11 the histogram behind it.

\(n\) is the index's k-mer count, which stands in for the residues in the database. --ka-k overrides \(K\).

In an hp alphabet a chance match rate near 0.40 is a polar protein against a membrane protein. They almost never match by luck, so the same \(S\) is rare between them. 0.64 is two membrane proteins, which match on most positions for free. That line is nearly flat because 0.64 sits just under the boundary \(2/3\), where \(\lambda_{\text{pair}}\) reaches 0.

1. Sample 200 sequences from the database.
2. Search them as they are, and search the same sequences shuffled so that every neighbouring pair of residues stays as common as it was (the reference).
3. Bin every region from both searches by \(x = \lambda_{\text{pair}}\,S\) in half-nat bins and count.
4. Per bin, take the log ratio of the two counts, \(r(x)\).
5. Fit a straight line through the eight lowest bins that both searches fill with at least 30 regions. That is the baseline; \(\sigma\) is the scatter of those bins around it. A constant ratio is a difference in \(K\); a gentle tilt is sequence structure the reference lacks.
6. Walk up in \(x\). The first bin whose ratio sits more than \(2\sigma\) above the line is where the real queries have something the same sequences, shuffled, cannot have: relatives. The fit stops there.
7. Fit \(\ln c_{\text{real}}(x)\) on the eight highest bins that made it in: \(r_{\text{database}}\) is minus the slope, \(K\) comes from the intercept.

In \(E = K \cdot m \cdot n \cdot e^{-\lambda_{\text{pair}} \cdot r_{\text{database}} \cdot S}\) the product \(m \cdot n\) is the number of squares in the grid of query positions against target positions, and the formula treats every square as a fresh chance at a region. That overcounts: squares next to each other on one diagonal belong to the same excursion, and a region is counted once, at its peak. \(K\) is the fraction of those \(m \cdot n\) squares that are real, independent starting points, so it is well under 1. For plain ungapped scoring Karlin and Altschul (1990) give a formula for \(K\) from the score distribution. kmerseek cannot use it, because its regions are not plain excursions. A region exists only if it contains an exact 12-mer seed, the give-up margin \(X\) ends it, and walks on one diagonal are merged. Each of those removes chance regions, mostly low-scoring ones, and the formula knows nothing about them. So \(K\) is measured instead, the way BLAST measures its gapped constants by simulation: it is what is left of the intercept of Equation (7a) once the sizes and the bin width are divided out, Equation (7c):

The seed requirement is why kmerseek's values (0.0227 on SCOPe40, 0.0180 on Swiss-Prot, 0.0139 on UniRef50) are several times smaller than BLAST's ungapped \(K\): most squares of the grid can never start a region because no 12-mer matches there. \(K\) is fitted once per database, pooled over every pair, so it is a region-weighted average; the theory says \(K\) varies with composition too, but far less than \(\lambda\) does, which is why the composition goes into \(\lambda_{\text{pair}}\) and not into \(K\).

Two things pull the fitted slope under the true average. Relatives enter gradually, so the bins just below the stop already carry a few. And a mixture of queries with different \(r_{\text{database}}\) is shallower than their average at high \(x\), because the shallow ones dominate there. The equation at the top of the page uses the constants measured for the chosen database and does not move with these sliders.

The two figures below are not kmerseek's calibration output. They come from a small simulation inside this page that makes up region counts with the shapes described above, so you can change things the real run does not let you change and watch what each one does to the fitted slope. In the simulation each calibration query departs from independent positions by its own amount, drawn around the r_database you set with the spread you set, and has its own share of relatives. The pooled counts are the sum over queries, so the r_database the fit recovers is a region-weighted average over them. Relatives add a shallow tail that begins around 9 nats, the reference is the chance curve without relatives, and Poisson noise is added to every count. The switches change what the reference keeps, how many queries you spend, how much they differ, and how the stop rule is set. The readouts compare the fitted slope with the true average the simulation was given.

\(r_{\text{database}}\) is fitted on chance regions only: the bins below the stop rule, where the real queries and the same queries shuffled agree. Relatives never enter those bins, so \(\kappa\), which describes relatives, cannot tilt them; it can only decide where the fit stops. What tilts them is run structure inside the sequences. Real proteins have helices, strands and membrane stretches where one class repeats for several residues in a row. The Karlin–Altschul model assumes each position is drawn on its own.

Figure 9 puts one knob on that. \(\rho\) (rho) is the neighbour repeat rate: at each position, with probability \(\rho\) the sequence keeps the class of the position before, and otherwise draws a class from its composition. Both sequences of a pair get it, so the chance that two neighbours carry the same class is \(\rho + ((1 - \rho) \cdot \Pr(\text{match} \mid \text{unrelated}))\), the same shape the coin of section 1 gives a match inside a shared domain, with \(\rho\) in place of \(\kappa\). At \(\rho = 0\) the positions are independent and the fitted slope is close to 1, the Karlin–Altschul value. When neighbours repeat, a lucky match tends to be followed by another lucky match, chance regions run longer, the tail of the counts falls more slowly, and \(r_{\text{database}}\) drops below 1. The shuffled reference keeps the neighbour pairs, so it keeps \(\rho\) and the stop rule still fires only on relatives.

Drag \(\rho\) and the whole tail tilts; the purple curve on the right is that tilt. Drag \(\kappa\) and the chance bins do not move, but the pink line does. Relatives pile up at the score their shared domain can reach: at \(\kappa = 0.46\) that is just above the chance regions, so the stop rule fires early and the fit is cut off where the tail is still bending, which is why the coral dot sits under the curve. At \(\kappa = 0\) there is nothing to stop on, and at \(\kappa = 0.9\) the relatives climb past every chance region, so in both cases the fit runs to the end of the counts. This is the same effect the readout of Figure 7 calls relatives entering gradually. The alphabet and \(C\) from section 2 and \(X\) from section 1 also enter, because they set which chance regions exist; the sliders of section 2 do not, because this model draws both compositions at the Swiss-Prot average.