How many independent proteins are in a benchmark? Error correlation across sequence similarity and the cost to error bars

Abstract

When a model misclassifies one protein it tends to misclassify its relatives, so a benchmark of n proteins holds fewer than n independent observations. We measured how much fewer. On 9,000 Swiss-Prot proteins from twelve organisms we aligned every pair, then estimated the correlation of 0/1 prediction errors as a function of Smith-Waterman similarity. It is positive in every similarity class that passes a significance gate, 0.01 to 0.07 for the weakest, and it sums over so many weak pairs that the effective size of a random subset saturates near 1,250 to 4,650 proteins. In a simulation with known truth, naive 95 percent intervals covered 0.78 at n = 1,600, clustering-based intervals 0.79 to 0.83, and a similarity-weighted interval 0.95.

Motivation

A test set of 1,000 proteins and a classifier whose error rate is 27.5 percent. The textbook 95 percent interval for the accuracy is plus or minus 2.8 percentage points. The enzyme-or-not classifier we study below has that error rate on such a set, and our estimate of its design effect on a random 1,000-protein subset is 1.80, which widens the interval to plus or minus 3.7 points. At 3,000 proteins the textbook interval is plus or minus 1.6 points and ours is plus or minus 2.9. Reporting the first number is not dishonest, but it quietly assumes that 3,000 proteins are 3,000 independent draws.

They are not, and everyone who works with protein benchmarks knows it. Homology-aware splits exist to stop leakage between training and test sets. What has received less attention is the dependence inside the test set. If a family of forty related kinases sits in the benchmark, the model's performance on that family is one fact reported forty times. The standard remedy in other fields is to cluster the units and resample clusters. That requires a threshold, and for proteins there is no threshold at which clusters are both independent and non-degenerate, as we show.

The question we ask is practical. Given a benchmark and a model, how many independent proteins does the benchmark contain, and what should the error bar be?

Related work

Splitting sequence families so that training and test sequences are at most a given identity apart is well developed. Petti and Eddy give independent-set algorithms that split Pfam families under identity limits on train-test and test-test pairs, and they down-sample their test sets to limit the weight of large families [1]. Their subject is construction; as far as we could determine they do not analyse the dependence among test sequences that remains. Cluster-aware intervals are used in applied protein work, for example bootstrap resampling over homology clusters at 40 percent identity [2], and cluster-adjusted effective sample sizes appear in the evaluation of classifiers on nested text data [3].

The statistics we use come from elsewhere. The design effect of clustered samples and the effective sample size n divided by it are due to Kish [4]. Phylogenetic effective sample size defines the number of independent observations carried by a sample with a known correlation structure [5]. Leinster and Cobbold's similarity-sensitive diversity measures count effective species from a similarity matrix [6]. In spatial econometrics, Conley's estimator weights residual cross-products by a kernel of distance instead of using clusters [7]; with a uniform kernel and a cutoff it becomes the estimator we use below. Henikoff and Henikoff down-weight related sequences when building profiles [8].

What is new here, to our knowledge, and only this: (a) a measurement of the correlation of model errors between proteins as a function of alignment similarity, on real proteins and real classifiers, showing that it is positive down to the weakest significant hits; (b) a demonstration, against a known truth, that naive intervals and cluster-robust intervals at similarity thresholds both lose coverage as benchmarks grow, because the dependence is carried by weak pairs that no usable clustering threshold joins; (c) an effective-size estimator calibrated on that correlation function, and a calibration-free pairwise estimator in Conley's style with an alignment-significance gate as the cutoff, both restoring coverage; (d) a sensitivity analysis showing that the answer depends on the significance gate and cannot be transferred between prediction tasks.

Method

Proteins and alignments

We took the reviewed Swiss-Prot entries of twelve organisms from the UniProt REST API (release 2026_03, accessed 2026-10-05) [9], restricted to 60 to 700 standard residues, and drew 1,000 proteins per organism with a fixed seed, giving a pool of 12,000. The organisms are human, mouse, budding and fission yeast, Escherichia coli K-12, Bacillus subtilis, Arabidopsis, Drosophila, Caenorhabditis elegans, Mycobacterium tuberculosis, Dictyostelium and zebrafish. We aligned all 71,994,000 pairs with local Smith-Waterman alignment [10], BLOSUM62 [11], gap open 11 and gap extend 1, using the Opal SIMD implementation through its Python bindings, whose documentation cites [12]. This took 1,672 s on three threads.

We converted each raw score S for proteins of lengths m and n into an expectation E = K m n exp(-lambda S), with lambda = 0.267 and K = 0.041, the constants the NCBI BLAST source tabulates for this matrix and gap cost. These constants are for database searches and are only an approximation for single pairs, so we treat E as a ranking device, not a probability. We define the similarity of two proteins as zero if E exceeds a gate of 1e-3, and otherwise as the bit score of the alignment divided by the smaller of the two self-alignment bit scores, a number between 0 and 1. About 539,000 of the 71,994,000 pairs (0.75 percent) pass the gate. At similarity 0.1 or above there are 85,610 pairs, at 0.3 or above 6,942, at 0.5 or above 2,047.

Classifiers and errors

We drew a random training set of 3,000 proteins and used the other 9,000 as the evaluation universe U. Two labels came from UniProt annotations: whether an entry has a transmembrane feature (tm), and whether it has an EC number (enz). Two kinds of predictors per label: logistic regression on amino-acid and dipeptide composition (420 standardised features, C = 0.05; "lr"), and nearest-neighbour label transfer from the best-scoring training protein under our similarity ("nn", with a best hit available for 97.9 percent of U). This gives four error vectors e (1 for a misclassified protein): tm_lr (error rate 11.7 percent), enz_lr (27.5), enz_nn (19.8) and tm_nn (14.0). We also use the difference of two prediction errors as an outcome for paired comparisons. The nearest-neighbour predictors use the same alignment scores as the similarity, so their error correlation is not independent of the measurement and we report them as a stress case.

Error correlation as a function of similarity

For a set of proteins with errors e_i, centre them, r_i = e_i - mean(e), and let v be the mean of r squared. We bin similarity into eight classes (no significant hit; (0, 0.05]; (0.05, 0.1]; (0.1, 0.15]; (0.15, 0.2]; (0.2, 0.3]; (0.3, 0.5]; (0.5, 1]) and estimate for each class b the correlation rho_b as the mean of r_i r_j / v over pairs in the class. We then force rho to be non-decreasing in similarity by isotonic regression weighted by pair counts, and clip at zero.

Given rho, the design effect of a benchmark B of n proteins is the sum of rho(s_ij) over all ordered pairs in B including i = j, divided by n; the effective size is n over that; and an error bar is the naive standard deviation multiplied by the square root of the design effect. For random subsets drawn from a universe of N proteins the expected design effect is 1 + (n - 1) times rho-bar, where rho-bar is the average of rho over all pairs of the universe; the effective size then approaches 1 / rho-bar as n grows, which we call the ceiling.

Centring by the sample mean removes any correlation that is shared equally by all pairs, so rho is identified only relative to the unrelated pairs. Any common factor that shifts all proteins together is invisible here and is not counted. The ceilings are therefore upper bounds on the information in the benchmark.

Tests of the interval estimators

We evaluated five interval estimators on a benchmark: naive (variance of the mean error from the sample variance divided by n); cluster-robust, grouping the benchmark's proteins by single-linkage components at a similarity threshold tau (we tried 0.1, 0.3 and 0.5) and using the standard sandwich variance over components; effective-size (calibrated rho as above); a Conley-style kernel estimator with no calibration, whose variance is the sum of r_i r_j over pairs that pass the gate (and i = j), divided by n squared; and, for the real-data validation only, an "exact" variant that also corrects for the benchmark being a part of a finite universe.

Real replicate benchmarks. Single-linkage components of U at similarity 0.3 (or 0.5) were randomly split into a calibration half and a validation half. We estimated rho on the calibration half. We then divided the validation half into K groups by assigning whole components at random to the currently smallest group (K = 24, 12, 6, 3, giving mean group sizes 187, 375, 750 and 1,500). Every group is a benchmark. Each estimator gives a standard error for its group's mean error; we score z = (group mean minus mean over the validation half) divided by that standard error, by the standard deviation of z (1 if calibrated) and by the coverage of a 95 percent interval. We repeated this with 4 random calibration and validation splits and 40 random group assignments per setting (3 and 25 for the kernel comparison).

Known-truth simulation. We took the calibrated rho for enz_lr from the calibration half of one split, built the correlation matrix it implies for 4,500 proteins plus 4,500 calibration proteins, clipped its negative eigenvalues (687 of 9,000) and rescaled to unit diagonal, and drew 300 independent Gaussian fields on those proteins. For each field we estimated rho from the calibration proteins alone, then drew 20 random subsets of each size from the evaluation proteins and scored each estimator against the true mean of zero. A second simulation of 200 fields and 15 subsets per size compared the calibrated and kernel estimators. The simulated outcome is continuous and Gaussian, not a 0/1 error; what is real is the similarity structure.

Null experiment. We shuffled the error vector among proteins, which destroys any link between similarity and error, and repeated the real-data procedure.

Results

Errors are correlated at every significance level

Across the 9,000 proteins there are 297,313 unordered pairs that pass the gate, among 40.5 million. The estimated error correlation rises monotonically with similarity for all four predictors and is not zero at the bottom.

The first row is not a result but a check on the centring: pairs without a significant hit show correlation within 0.001 of zero. In the bottom classes correlations of 0.01 to 0.12 look small, but there are 248,000 such pairs, against about 9,500 pairs with similarity above 0.2. The product dominates: for enz_lr, pairs with similarity at or below 0.1 contribute 0.073 x 130,656 + 0.123 x 117,625, which is about 24,000 units of correlation, against about 8,500 for all pairs above 0.1.

Why clustering does not rescue the interval

Single-linkage components on the pool of 12,000 are degenerate at any threshold that would capture the weak pairs.

At 0.1, a single component holds 91 percent of the pool, and so a cluster-robust interval at that threshold has one cluster that is the whole benchmark. At 0.3 or 0.5 the components are small and honest, but they join only the 6,942 and 2,047 strongest pairs, not the quarter-million weak pairs where most of the correlation sits.

How many independent proteins?

Applying the calibrated rho to random subsets of the 9,000-protein universe gives design effects and effective sizes that depend strongly on the task.

The cells give the design effect. Equivalent effective sizes for enz_lr are 93, 242, 555, 882 and 1,096 proteins at the five sizes: the full 9,000-protein universe behaves like 1,096 independent proteins, and no larger benchmark of this kind can do better than about 1,250 for this task.

These figures describe a universe: 12 organisms, 750 evaluation proteins each, sampled uniformly from reviewed entries. Adding proteins from more families would lower rho-bar and raise the ceiling; a benchmark that over-samples a few large families would do the opposite. The number is a property of the benchmark and the model, not of proteins in general.

Stability. Recomputing rho-bar on 24 random halves of the proteins gave a mean within 1.2 percent of the full-universe value for enz_lr and a standard deviation of 6.7 percent (tm_lr: 17 percent, enz_nn: 7.8 percent, tm_nn: 19 percent). The 2.5 to 97.5 percent range of the half-sample ceiling was 1,123 to 1,451 for enz_lr, 2,070 to 3,750 for tm_lr, 1,410 to 1,840 for enz_nn and 3,370 to 6,400 for tm_nn. Halves chosen as whole 0.3-components gave means within 12 percent of the full value (ceilings 2,512, 1,185, 1,513 and 4,143 for tm_lr, enz_lr, enz_nn and tm_nn) with spreads between 13 percent smaller and 68 percent larger.

Transfer. The calibration does not carry over between tasks. If we applied the rho curve of one task to another, the resulting rho-bar was off by a factor between 0.27 and 3.7 (for example, enz_lr's curve applied to tm_nn: 3.7 times the true value; tm_nn's curve applied to enz_lr: 0.27). Between two predictors of the same label the factor was 0.6 to 1.7. A default curve is not safe; the curve should come from the model under evaluation.

Sensitivity to the gate. The ceiling is a function of which pairs are counted as related.

Between 1e-5 and 1e-1 the ceiling moves by a factor of six for enz_lr and five for tm_lr. At 10, 99.9 percent of pairs pass, the centred sum over passing pairs is nearly the full sum, which is zero by construction, and the estimate is meaningless. If E behaved as a p-value for an unrelated pair, about a tenth of a percent of unrelated pairs would pass a gate of 1e-3 by chance, roughly a seventh of the 297,313 pairs that pass in U; we did not verify that calibration. Replacing the minimum of the self-scores by their geometric mean in the normalisation changed no ceiling, because the ceiling depends only on which pairs pass the gate.

Coverage when the truth is known

In the simulation the true superpopulation mean is zero, and we know the true correlation structure. Rows give the coverage of nominal 95 percent intervals (and, in brackets, the standard deviation of z) over 300 fields times 20 subsets.

Naive intervals lose coverage steadily as the benchmark grows, and the cluster-robust intervals recover only a few points of it: 0.81, 0.83 and 0.79 at n = 1,600 against 0.78. The calibrated effective-size interval stays at 0.95. Treating each of the 300 fields as an independent trial, the binomial standard error of a coverage value is about 0.013, so differences of 0.01 between the effective-size column and the one with the true correlation are within noise, but the 0.18 gap to naive at n = 1,600 is not.

The estimated rho from a single field was biased low: the mean estimated curve was 0.068, 0.093, 0.123, 0.170, 0.270, 0.325 and 0.504 in the seven significant-hit classes against true values 0.089, 0.121, 0.144, 0.193, 0.307, 0.342 and 0.520, so 3 to 24 percent too small, and the estimated design effect was 6 to 12 percent below the true one at the three largest sizes. The intervals still reached nominal coverage. We did not investigate why; the finite-sample variance estimate and the bias in rho may be offsetting each other.

Calibration-free kernel estimator

The second simulation compared the calibrated estimator with the kernel estimator that needs no calibration set: sum the products of centred errors over pairs that pass the gate.

The kernel estimator's variance was never negative in 3,000 subsets per size, and its coverage was within 0.01 of the calibrated estimator's at every size, a little lower at 200 (0.937 against 0.946) and 1,600 (0.935 against 0.943). For someone who has one benchmark and one model, the kernel estimator is the practical choice; the calibration is useful for prediction, where one wants the curve before running the next model, and for the ceiling.

Real replicate benchmarks

On real proteins we can only construct replicates that are separated at some threshold, and weak pairs necessarily connect them. At threshold 0.3, groups of 187 and 375 proteins gave these results for enz_lr (the other combinations behave alike; the kernel column comes from a separate run with 3 splits and 25 assignments, so its naive baseline differs slightly, 1.14 and 0.91 at 187 and 1.10 and 0.92 at 375):

For the four label-and-predictor combinations at group size 187, the naive standard deviation of z was 1.12 to 1.18 and coverage 0.91 to 0.92; the exact estimator gave 0.99 to 1.00 and 0.95 for all four, and the effective-size interval 0.94 to 0.98 and 0.96. With whole groups of 24 per partition and four calibration splits, we have roughly 100 independent groups per cell, so the binomial standard error of a coverage value is about 0.026: the gap between naive and exact at n = 187 is about 1.5 standard errors per cell, suggestive across four cells but not conclusive on its own.

At larger group sizes (750 and 1,500) every estimator over-covers, including naive, and the effective-size interval most. This is the expected behaviour of this design. A group of 1,500 is a third of the validation half, and its deviation from the mean over that half is shrunk by construction, whereas the effective-size interval targets deviation from a superpopulation. The "exact" variant, which corrects for this, had coverage 0.98 and 0.99 at 750 and 1,500 for enz_lr and standard deviation of z 0.87 and 0.81, i.e. still conservative. We cannot validate effective-size intervals on real data at large n, which is why the known-truth simulation matters. At threshold 0.5 the groups are less separated, the observed design effect at group sizes 187 and 375 was only 1.04 to 1.16, and naive intervals covered 0.93 to 0.96.

Running the real-data procedure with gates 1e-5 and 1e-1 gave the same picture for the exact estimator (coverage 0.95 to 0.96 at 187 and 375) and a more conservative effective-size interval at 1e-1.

Paired comparisons and the null

The same machinery applies to the per-protein difference between two models' errors. For enz_lr against enz_nn the observed design effect between groups was 1.27 at 187 and 1.34 at 375, with naive coverage 0.93 and 0.94 and effective-size coverage 0.96 and 0.97. For two exchangeable models (logistic regressions trained on random halves of the training set, so that any difference is training noise) the observed design effect was 1.29 and 1.37; even these differences are correlated across relatives, in the sense that the way two equally good models disagree has family structure.

The shuffled-errors null behaved as it should. The estimated rho was within 0.02 of zero in every bin, the estimated design effect was 1.01, the observed one 1.01 and 1.00, and naive, effective-size and cluster-robust intervals all covered 0.95 to 0.96. The correction does not inflate error bars when there is nothing to correct.

Falsification and limits

We tried to break the idea in four ways, and some things did not survive.

The strongest claim of the real-data design failed to materialise. We hoped that cluster-disjoint replicate groups would give an independent test of the effective-size interval. They do not, for a reason that is itself a result: groups separated at a threshold remain linked by weak pairs, which are the pairs that carry the correlation. So at group sizes of 750 and above the design cannot discriminate between estimators, and the known-truth simulation carries the argument for large benchmarks. That simulation has its own weaknesses. The outcome is Gaussian and continuous, the correlation matrix required clipping 687 negative eigenvalues (so the realised covariance is slightly smaller than the one we assumed, which is why the oracle intervals are mildly conservative), and the correlation function is the one we estimated, so it cannot fail to resemble the data in its shape.

The gate is a hidden degree of freedom. We fixed 1e-3 before any evaluation, but the ceiling moved sixfold across reasonable gates, and the one choice we can defend, that unrelated pairs should not enter, is violated at 1e-3 by an estimated one seventh of the passing pairs. The ceilings should be read as indicative orders of magnitude, with the 1,250 to 4,650 range for the gate we used. The E-values rely on search-database constants for pairwise comparisons.

Centring hides any correlation that all proteins share. A benchmark drawn from one lab's annotation pipeline, or all from one clade, can carry shared errors that this approach cannot see. The ceilings are upper bounds.

The experiments are on one universe (twelve organisms, Swiss-Prot annotations, composition-based and nearest-neighbour predictors, two labels). We do not know how the curves look for protein language model embeddings, structure prediction, or labels measured by experiment. The nearest-neighbour predictor shares its alignment with the similarity, so its numbers are a stress case. The mean 0/1 error is one metric; we did not look at AUC or at regression targets. Coverage differences between 0.91 and 0.95 on real replicates are within 1.5 standard errors per cell.

What this means in practice

For a benchmark of a few hundred proteins, the textbook interval is close to right. In our universe, at 300 proteins, the design effect was 1.06 to 1.24 depending on the task, which widens an interval by 3 to 11 percent. The error bar becomes materially wrong when the benchmark has thousands of proteins, because the design effect grows roughly linearly with n; beyond about 1,000 to 3,000 proteins in this universe, adding more proteins narrows the interval very little.

With one model and one benchmark, use the kernel estimator: align the benchmark against itself, gate at E of 1e-3 (report the gate), and sum the centred error products over pairs that pass. In a known-truth test this gave 0.935 to 0.947 coverage at 400 to 1,600 proteins where the naive interval gave 0.76 to 0.90. Report the number of proteins and the effective size next to the accuracy. Do not use a clustering threshold as a stand-in without checking that it does not collapse into one giant component, and do not borrow a correlation curve from another task.

For comparisons between two models, apply the same estimator to the per-protein difference of errors; do not assume that the difference is independent because the models are similar.

Reproducibility appendix

Seeds. Pool: 1,000 proteins per organism after a shuffle with seed 20261005. Training and evaluation split: seed 20261006 (3,000 and 9,000). Calibration and validation splits: seeds 1000 to 1003; random group assignments 10000 x split + index. Known-truth simulation: split seed 4242, random generators 99 (300 fields) and 2024 (200 fields). Half-sample stability: seed 5. Times: alignment 1,672 s on three threads (x86_64 SIMD under Rosetta); real-data validation 800 to 1,500 s per setting under heavy machine load; simulations 2,622 s and 758 s.

import numpy as np
from sklearn.isotonic import IsotonicRegression

LAM, KK = 0.267, 0.041
EDGES = np.array([0, 1e-9, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 1.01])

def similarity(S, selfs, lengths, gate=1e-3):
    bits = lambda x: (LAM * x - np.log(KK)) / np.log(2)
    E = KK * np.outer(lengths, lengths) * np.exp(-LAM * S)
    s = np.clip(bits(S) / np.minimum.outer(bits(selfs), bits(selfs)), 0, 1)
    s[E > gate] = 0.0
    np.fill_diagonal(s, 1.0)
    return s

def bins_of(s):
    b = (np.digitize(s, EDGES) - 1).astype(np.int8)
    np.fill_diagonal(b, -1)
    return b

def calibrate(B, y):
    r = y - y.mean()
    v = (r ** 2).mean()
    rr = np.outer(r, r)
    nb = len(EDGES) - 1
    cnt = np.array([(B == k).sum() for k in range(nb)], dtype=float)
    num = np.array([(rr * (B == k)).sum() for k in range(nb)])
    rho = num / np.maximum(cnt, 1) / v
    iso = IsotonicRegression(increasing=True, y_min=0.0, y_max=1.0)
    lut = iso.fit(np.arange(nb), rho, sample_weight=np.maximum(cnt, 1)).predict(np.arange(nb))
    lut[0] = max(0.0, min(lut[0], lut[1]))
    return lut

def design_effect(B, lut):
    R = lut[np.where(B < 0, 0, B)]
    np.fill_diagonal(R, 1.0)
    return R.sum() / len(B)

def kernel_se(e, s):
    r = e - e.mean()
    W = s > 0
    np.fill_diagonal(W, True)
    return np.sqrt(max((np.outer(r, r) * W).sum() / len(e) ** 2, e.var(ddof=1) / len(e)))

def effective_se(e, B, lut):
    return np.sqrt(e.var(ddof=1) * design_effect(B, lut) / len(e))

The Smith-Waterman scores for 12,000 proteins were computed as blocks of 500 sequences with pyopal.align(query, database, scoring_matrix='BLOSUM62', gap_open=11, gap_extend=1, algorithm='sw', mode='score').