The performance of genetic-constraint metrics varies significantly across the human noncoding genome

A longstanding goal in human genetics is to prioritize noncoding loci that, when disrupted, lead to developmental disorders and other Mendelian traits. In pursuit of this goal, multiple metrics have been developed to distinguish neutrally evolving sequences from those subjected to purifying selection. These metrics are commonly evaluated genome-wide, e.g., by computing a precision-recall curve on windows tiling the entire noncoding genome. Here, we identify parts of the noncoding genome where these metrics significantly underperform relative to their genome-wide performance due to “bias” in the underlying models of neutral genetic variation and/or a low “signal-to-noise ratio” in the genetic data. The most extreme effects are found for Gnocchi (Chen et al. 2024), the performance of which declines as GC content increases. We suggest annotating constraint scores of noncoding genomic intervals with robust measures of the bias of the corresponding model, allowing users to gauge confidence in those scores.

Computational approaches designed to predict noncoding variant pathogenicity must be improved if they are to have a meaningful impact upon clinical research (Liu et al. 2019). Many of these approaches are based on a simple, but powerful, idea: over the hundreds of millions of years separating species, a single base is likely to have mutated, so the only plausible way to explain why a genomic site is identical between a human and a fish, say, is that whatever mutations occurred there were harmful and eliminated by selection. This deep-time filter is so rigorous that it can detect even weakly constrained elements. Methods like phyloP (Pollard et al. 2010) and GERP (Davydov et al. 2010; Cooper et al. 2005) formalize this idea in a number of ways. (1) They go beyond the vague assertion that a base is "likely to have mutated" by quantifying the expected (neutral) rate of nucleotide substitution between species. (2) They extend the analysis from the comparison of two species to the entire tree of life by accounting for the fact that conservation between, e.g., human and chimp (closely related) is weaker evidence of selection than conservation between human and fish (distantly related). (3) They go beyond a simple binary classification of conservation to answer the question "Is the substitution pattern we see at this site extremely unlikely to have occurred under the neutral model?” Selection signals derived from sequence conservation are extraordinarily strong and unambiguous but have an Achilles heel: they cannot find genomic elements that are newly 1 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint constrained in humans, limiting their ability to correctly prioritize noncoding variants. If an element (like an enhancer) becomes functional in the human lineage after the split from chimpanzees, it will not be conserved in chimps, mice, or fish. In scenarios such as these, conservation methods are underpowered (though they can capture elements that have gained function in closely related species, e.g., primates (Boffelli et al. 2003)). The blind spot of conservation-based methods is not small. Rands et al. estimated that only about 25% of the elements under purifying selection in humans are also constrained in mouse (Rands et al. 2014), while, more recently, Huber et al. estimated that most of the noncoding sequence that is functional in human has experienced changes in selection coefficients throughout mammalian evolution (Huber et al. 2020). This process, called “turnover”, had long been suspected from the observation that, as the divergence between mammalian species increases, the predicted amount of pairwise shared functional sequence decreases dramatically (Smith et al. 2004; Meader et al. 2010). Furthermore, most transcription-factor binding events are species-specific, indicative of gain and loss of binding along each lineage (reviewed in (Ponting and Hardison 2011)). The only way to discover human-specific genetic constraints in the noncoding genome are population-based methods. Pleasingly, early deployment of these methods discovered that regulatory elements located near neural genes gained function in the human lineage (Ward and Kellis 2012; Schrider and Kern 2015), and we expect these methods to become increasingly sensitive for the discovery of gain- and loss-of-function discovery as more human genomes are sequenced. These approaches start with the construction of a model of neutral genetic diversity (e.g., number of SNVs in each of a set of 1kb genomic windows) observed in a healthy human cohort, using features that are known to co-vary with SNV count. Of these, the most well-known feature is the local sequence composition of a given site, whose probability of being polymorphic can increase >10-fold if the sequence is a CpG dinucleotide (Sved and Bird 1990; Bulmer 1986; Schraiber et al. 2024). With a model of this type in hand, one may identify putative signals of negative selection by searching for anomalous regions in the genome where the observed diversity is significantly lower than the neutral model predicts (Chen et al. 2024; Halldorsson et al. 2022; di Iulio et al. 2018). Alternatively, one may construct a model of all genetic diversity, including observed diversity impacted by negative selection. In this approach, the model includes a selection parameter, the value of which is estimated for each genomic window (Dukler et al. 2022). The two approaches are intimately related: as Dukler et al. show, the maximum-likelihood estimate of the selection coefficient in the second approach is the ratio of the observed diversity and that expected under neutrality, the comparison of which is the basis of the first approach. The success of either approach depends upon how well the underlying model of neutrality tracks the observed SNV density across the noncoding genome: if one fails to or inadequately models the neutral mechanisms that, like purifying selection, deplete SNV density, the result is a “biased” model whose predictions of genetic constraint will be flawed. Even if a neutral model can perfectly recapitulate the expected distribution of SNV counts in a locus, one's ability to predict that the locus is under purifying selection is still limited by how large the “signal” (the depletion of SNV counts due to purifying selection) is relative to “noise” (unavoidable statistical fluctuations in observed SNV counts). Motivated by these concerns, we investigate how existing 2 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint models’ ability to predict genetic constraint is impacted by variations in model bias and signal-to-noise ratio.

Consider a recent model of neutral genetic diversity in humans (Chen et al. 2024), which predicts the distribution of the number of SNVs in a given genomic window, aggregated over 76,156 humans, as a function of local sequence context, nucleotide-level DNA methylation, and a comprehensive set of regional genomic features, including GC content and recombination rates. This model yields a score called Gnocchi, denoted below, which evaluates how typical or unusual an observed SNV count is in the context of the window containing the SNVs, where the residual, , is defined by , and is the observed SNV count, and and are the mean and standard deviation, respectively, of the predicted distribution of SNV counts for that window. A value of Gnocchi near zero indicates that the observed SNV count is close to what one would expect from the neutral model, . Gnocchi is a standardized residual: the residual, , is reported relative to the expected background noise, . Accordingly, indicates the SNV count is not only smaller than , but also smaller than expected statistical fluctuations below , raising the possibility (but not guaranteeing) that negative selection is operative in that window. Standardized residuals provide information about how the underlying model of neutral diversity is biased. To illustrate how, we developed a simulation consisting of the following sequence of steps. (i) We simulated mutation by random generation of SNVs in a set of mock genomic windows as a function of a mock genomic feature (an example of a real “genomic feature” is “nucleotide content”) distributed normally over windows (Figure 1A). The “generative process” is characterized by a “true rate” at which mutations are generated as a function of the genomic feature (black line in Figure 1B-D), and sampling from a Poisson distribution with that rate (grey dots in Figure 1B-D). (ii) We simulated the development of the model of neutral diversity underlying Gnocchi by estimating the parameters of various models that map the feature value of each window to its SNV count. Like the generative process, the estimated model has its own rate, which we call the “learned rate” (orange line in Figure 1B-D) (iii) We simulated Gnocchi by combining the results of (i) and (ii) to compute a standardized residual (grey dots in Figure 1E-G). The difference between the estimated model and the generative process is quantified by “global model bias”, defined as the squared difference between the learned rate and the true rate, averaged over all windows (Bishop 2006). We started by using the generated SNVs to estimate a model of the same functional form as the generative process (the “Quadratic model”; Figure 1B). Consequently, after training, the estimated model recapitulated the simulation's generative 3 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint process, resulting in zero global model bias (Figure 1B). Corresponding to this lack of bias, the distribution of the standardized residual is randomly scattered around zero independent of genomic feature, indicating that observed and predicted SNV counts agree, on average (Figure 1E). If the estimated model has a simpler functional form than the generative process, such as when the estimated model is linear but the generative process is quadratic, then the global bias becomes non-zero (Figure 1C). Associated with the global bias is a “feature-specific bias”, defined as the mean of the residuals as a function of feature value (orange line in Figure 1F). The feature-specific bias shifts away from zero, particularly in the tails of the distribution of feature values, indicating that the model systematically over- or under-estimates observed genetic diversity. When the model is an even poorer approximation of the generative process, the global bias is even larger (Figure 1D), and the feature-specific bias deviates significantly from zero even in the center of the feature distribution (Figure 1G). Figure 1 Biased models over- or under-estimate SNV counts, on average, in the tails of the feature distribution. (Top row) Given a mock set of genomic intervals, each with a genomic feature (distributed as shown in A) and a corresponding SNV count (y-axis in B-D), we used Poisson regression to fit three different models mapping the feature to the SNV count (B-D). The quadratic model fits the data perfectly as SNV counts were generated by sampling from a Poisson distribution with a rate parameter that is a quadratic function of the genomic feature (note the concordance of the black line and orange line in B). The linear and constant models of the rate parameter are progressively less complex than the true rate parameter. In these cases, minimization of the sum of negative log likelihoods over all the windows during training forces the learned models to pass through the observed data where they are dense at the expense of deviating systematically from them where they are sparse: compare the grey data points and the orange line in panels C and D. (Bottom row) The degree to which the models over- or under-estimate SNV count can be quantified by the standardized residuals, shown in grey. The average of these data points (orange line) show the feature-dependent biases that emerge when learned models are incomplete. The simulation can be found at the following URL: http://github.com/quinlan-lab/constraint-tools/blob/main/papers/neutral_models_are_biased/9.regression/simulation.4. 2.1.interplay-of-SNR-and-bias.ipynb. 4 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Motivated by these results, we sought to visualize the bias of the Chen model by plotting the distribution of Gnocchi on a set of 693,270 putatively neutral 1kb intervals (Methods) as a function of three features (Figure 2A-C): a measure of nucleotide composition (GC content); a measure of background selection (BGS) computed by (Murphy et al. 2023); and a measure of GC-biased gene conversion (gBGC) computed by (Glémin et al. 2015). We chose GC content as single-nucleotide substitution probabilities have long been known to decline as GC content increases (Fryxell and Moon 2005; Arndt et al. 2005; Mugal and Ellegren 2011). Additionally, there is evidence that GC content at kilobase scale affects the probability of a CpG site mutating, as the deamination process requires local DNA melting, which in turn requires more thermodynamic energy for bonds between G and C nucleotides (Elango et al. 2008). Finally, even the GC content of the local sequence context of a CpG site may correlate with its mutation rate (Chandra and Gao 2025). Though BGS and gBGC metrics are not explicitly included in the Chen model, they are neutral mechanisms well-known to affect diversity in the human genome (McVicker et al. 2009; Pouyet et al. 2018; Murphy et al. 2023), and have been modeled when inferring natural selection in the context of population genetics, e.g., (Telis et al. 2020). BGS, as a process, reduces variation in regions linked to sequences constrained by purifying selection (for a review see (Charlesworth 2013)). GC-biased gene conversion (reviewed in (Duret and Galtier 2009)) describes how mismatch repair of heteroduplexes formed during meiosis is biased towards GC alleles, and thus against AT alleles, altering their probability of fixation. Therefore allele frequency is also altered in a manner that mimics natural selection though is neutral (Duret and Galtier 2009). As such, like GC content, both BGS and gBGC are neutral processes that lead to reduced genetic diversity which could masquerade as selection to models that fail to account for them. 5 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Figure 2 Visualizing the bias of the Chen model of neutral genetic diversity. (A-C) The standardized residuals of the Chen model (Gnocchi) deviate from zero, on average, when the genomic windows they are assessed on exhibit extreme values of GC content, background selection (BGS) or GC-biased gene conversion (gBGC). Since we report the (standardized) rank of Gnocchi on the y-axes, the marginal distribution of Gnocchi (the distribution of Gnocchi over all windows, irrespective of feature value) is uniform with an average value of 0.5 (horizontal black lines). The dark-grey lines are computed by conditioning on each feature bin before computing the mean Gnocchi rank. Note that the slope of the dark-grey lines denote the marginal (not independent) effects of each feature on bias. For example, a position on the dark-grey line in (B) equals the value of Gnocchi conditioned on a specific value of BGS, but integrated over all possible values of GC content and gBGC. The dark-grey lines deviate from the expected value of 0.5, particularly in the tails of the feature distributions (where the heat map is dark). However the dark-grey lines lie at a height of approximately 0.5 when conditioned upon the mean of the feature distributions (vertical black lines). For (B), this implies that the Chen model captures the expected reduction in SNV counts arising from the weak levels of BGS present in the majority of windows, even though the model does not explicitly include BGS, presumably because it covaries with features that the model does include. The co-variation of features with one another allows the linear regression (light-grey line) to track the observed nonlinear behavior (dark-grey line) in (B). (D) Even in the vicinity of mean GC content, the standardized residuals of the Chen model (Gnocchi) exhibit greater variance than expected (from a standard normal distribution), indicating the existence of neutral diversity not captured by the model. GC content was chosen to be within 0.1 * std(GC_content) of the mean GC content, where “std” refers to “standard deviation”. The curve marked “Expected” was constructed by scaling a standard normal distribution according to the bin width of the histogram and the number of intervals. An explanation of why a standard normal distribution is expected using a simplification of the Chen model can be found here: https://github.com/quinlan-lab/constraint-tools/blob/main/papers/neutral_models_are_biased/10.residuals-are-over-di spersed.ipynb. 6 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint We observe bias (represented by the deviation of the dark-grey line from y = 0.5 in Figure 2A-C) in the Chen model in rare windows that harbor extreme levels of background selection, or GC-biased gene conversion, but particularly those windows with extreme levels of GC content. Note, however, that the lack of bias in the Chen model at intermediate feature levels does not imply that the model cannot be improved there. In fact, we observe that the variance of Gnocchi in that vicinity is larger than what would be expected (Figure 2D), implying that there remains variation in SNV count that the Chen model is unable to capture even at genomic windows whose feature values are observed frequently across the genome. We next sought to measure how each feature independently contributed to the Chen model’s bias. The slopes of the dark-grey lines in Figure 2A-C do not measure these contributions. Instead, these slopes represent not only the direct effect of any given feature (e.g, GC content) on model bias, but also the indirect effect of the complementary features (BGS metric and gBGC) due to correlation among these features. To capture the independent effect of each feature on model bias, we fit a linear model of all features to Gnocchi (Methods). The estimated parameters of the linear model represent the direct effects of the corresponding feature on the bias of the Chen model (Table 1). We find, for example, that the direct effect of the BGS metric on Gnocchi is negative, implying that stronger background selection (smaller values of the BGS metric) is associated with larger values of Gnocchi. This accords with our naive expectation that a model of neutral diversity that does not explicitly account for BGS (as is the case with the Chen model) ought to overestimate diversity where background selection is strong (assuming that BGS does not significantly co-vary with features that the Chen model does include). Interestingly, though a positive value of Gnocchi erroneously predicts constraint in the neutral windows subject to strong BGS, it correctly predicts constraint nearby, by definition of BGS. Since we standardized all features prior to fitting our linear model to Gnocchi, we may compare the values of the corresponding regression coefficients, revealing that GC content has the largest effect (in magnitude) on model bias: GC content explains almost as much variance in Gnocchi as explained by the three features jointly (compare for one feature vs three features, respectively, in Table 2). Table 1. The direct effect of each feature on the bias of the neutral models underlying each constraint metric. The regression coefficients may be compared not only within a given model, but also among models, since the features are standardized prior to fitting, and the target variables (ranks) in the fitting process are uniformly distributed between 0 and 1. 7 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Table 2. The proportion of variance in each constraint score that is explained, , by all three features considered in this study, and by the single feature with the largest direct effect on bias. To assess the generality of these results, we repeated our analysis for three other constraint metrics: (Dukler et al. 2022), Depletion Rank (Halldorsson et al. 2022) and CDTS (di Iulio et λ𝑠 al. 2018) (Methods). The results, presented in Tables 1 and 2, and Supporting Figure 1, are similar to those obtained for the Chen model in the sense that all models develop biases as one moves into the tails of the feature distributions, though the details differ. For example, whereas high GC content biases Gnocchi above zero (rank > 0.5), it biases below zero (rank < 0.5) λ𝑠 (Supporting Figure 1). These high-GC-content intervals may be enriched for promoter regions, where the authors of the metric observed a similar bias, which they attributed to λ𝑠 misspecification of the underlying neutral model there (Dukler et al. 2022). Another difference 8 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint between the Chen model and the others is that BGS, not GC content, explains most of the bias observed in Depletion Rank and CDTS (e.g., Table 2). It is unclear why, but part of the answer may lie in the fact that the Chen model includes GC content as a factor, but the models underlying Depletion Rank and CDTS do not. The procedure we adopted to identify the feature that drives model bias is closely related to Gradient Boosting (Friedman 2001), a technique for improving a model by regressing its residuals onto its features, and then adding the resulting “residual model” to the predictions of the original model. Viewed this way, the values reported in Table 2 measure how much an improved modeling of the corresponding features (e.g., via Gradient Boosting) is expected to improve the underlying model of neutral diversity. Comparing the values among features and models leads to two interpretations. First, the largest improvement to the Chen model is expected to come from more complex representations of sequence composition (for which GC content is a surrogate), whereas the largest improvements to the models underlying Depletion Rank and CDTS are expected to come from better representations of BGS. Second, of all the models, the Chen model stands to benefit most from an improved modeling of the three features considered in this study. The purpose of constraint metrics, such as Gnocchi, is to discover “genetic constraint”, i.e., to distinguish genomic windows that are neutral from those that are under purifying selection. To understand how well this approach is likely to work under various circumstances, we returned to simulation. The standard approach to modeling purifying selection is to use forward-time population genetic simulations, e.g., via stdpopsim (Adrion et al. 2020; Lauterbur et al. 2023; Gower et al. 2025). Though stdpopsim incorporates a host of neutral processes, each of them is expected to contribute to “noise” (unavoidable statistical fluctuations in observed SNV counts in genomic windows), obscuring the patterns we wish to illustrate with our simulations. To minimize this noise, and therefore maximize our ability to recover “signal” (the depletion of SNV counts due solely to purifying selection), we decided instead to return to our simple simulation of neutral diversity (Figure 1A-D), which is stripped of most of the neutral processes included in stdpopsim. We extended the simple simulation to model purifying selection by simply reducing the SNV counts of a random subset of our simulated windows, labeling them “constrained”, and labeling the rest “neutral”. We then built a constraint predictor by thresholding residuals corresponding to each of the models of neutrality depicted in Figure 1B-D: windows with residuals greater than the threshold were predicted to be constrained, whereas those with residuals lower than the threshold were predicted to be neutral. Finally, we assessed the performance of those models by comparing predicted labels with actual labels. Even if a model is unbiased, its performance is limited by the signal-to-noise ratio (SNR): the ratio of the “signal” we are trying to recover (the depletion of SNV counts due to purifying selection) and the unwanted background “noise” we must contend with (unavoidable statistical fluctuations in observed SNV counts). To illustrate this problem and the effect of variations in SNR along the genome, we conditioned on windows with three different values of the genomic feature (Figure 3A-C), and in each case computed the SNR (Figure 3D-F; Methods). Regardless of model bias, performance declines (Figure 3G-I) as SNR declines (reflected in an 9 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint increasing overlap of the residual distributions of the neutral and constrained windows; Figure 3D-F). Performance is impacted not only by low SNR, but also by high model bias. To deconvolve these effects, we introduced, for each of the three cases in Figure 3A-C, enough variability in feature values to expose the impact of bias on the performance of the three models, yet not enough to alter the SNR in each case (compare narrow with wide red strips in Figure 3A-C). With this change, not only does performance decline due to declining SNR as found previously (notice the similarity of the green lines in Figure 3J-L with the corresponding green lines in Figure 3G-I), but also because of increasing bias for any given SNR (for each panel of Figure 3J-L, compare the three colored lines). Figure 3 SNR and model bias are both expected to significantly affect the ability to call constraint. The narrow red strips in panels A-C define the subsets of simulated genomic windows that were used to compute panels D-I. Each subset contains neutral (orange) and constrained (blue) windows with residuals distributed as shown in panels D-F (with respect to the unbiased quadratic model). Also shown in those panels are the corresponding values of SNR, computed from the the standard deviation of the neutral residuals and the depletion of the simulated SNV counts in the constrained windows (Methods). The SNV-count depletion ( ) and class-imbalance (relative numbers of “constrained” and “neutral” windows) were chosen to be representative of real data (Supporting Figure 2). For each window, and each model, we predict that the window is constrained if the residual exceeds a threshold; by varying the threshold we sweep out a PR curve for each model, showing the independent effect of SNR on model performance (panels G-I). To see the independent effect of model bias, compare the PR curves for each model in panels J-L, 10 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint which correspond to the wide red strips in panels A-C, respectively. “Local bias” is defined similarly to the definition of “Global bias” in Figure 1, but this time using only those genomic windows that comprise the corresponding red strips in panels A-C. The distributions of the constrained windows in panels D-F peak at , which measures the depletion of SNV count that we implemented in the simulation, which can be found here. Supporting Figure 3 shows similar results when the class distributions are computed using standardized residuals, instead of raw residuals, as shown in panels D-F. Supporting Figure 4 shows similar results when SNV counts are depleted by an amount proportional to the feature value (i.e., is now feature-dependent), except that now SNR increases (instead of decreasing) as one moves into the right tail of the feature distribution. Given that model bias and SNR are expected to jointly affect the ability to predict constraint in a way that is impossible to predict without knowledge of that which we are seeking (genetic constraint), we set out to empirically measure performance of Gnocchi as a function of GC content, BGS and gBGC. We designed a simple classifier that predicts that a window is constrained if Gnocchi exceeds a threshold. Since the main purpose of, and main challenge facing, scores such as Gnocchi is to identify constraint in the noncoding genome, an example of which may be a conserved enhancer, we constructed a set of noncoding windows and labeled those that overlap a GeneHancer enhancer (Fishilevich et al. 2017) “constrained”, and those that don’t “neutral” (Methods). Measured against this truth set, Gnocchi clearly does better than random guessing (compare black line and dashed line in Figure 4A), but its performance in GC-rich windows is significantly reduced (red lines); conversely, in GC-poor windows, performance is better than expected (blue lines). Since Gnocchi is biased at low GC content (Figure 2A), and in light of the results in Figure 3, we conjecture that the increased performance at predicting constraint in GC-poor regions is due to an increased SNR. Similar variability in performance was observed as a function of BGS and gBGC (Supporting Figure 5A, B). These results illustrate how model bias and/or SNR jointly generate significant variation in ability to discover constraint under the Chen model, depending upon the nature of the windows being interrogated. 11 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Figure 4 The ability of models of neutral diversity to predict constraint varies with the nature of the genomic windows. (A) The performance of the Chen model at predicting constraint (on a lax truth set defined by whether windows overlap GeneHancer enhancers) can be made to vary significantly by changing the nature of the windows upon which performance is assessed. We varied the Gnocchi threshold to trace out the precision-recall curves. As a benchmark, the performance of a classifier that randomly guesses whether a window is constrained or not has a precision pegged to the fraction of windows that are “positive”, i.e., labeled “constrained” (dashed line; see Methods). To compare performance across GC-content bins, which vary in the proportion of windows labeled “constrained”, we downsampled the bins so that each had an equal fraction of “constrained” windows. (B) Performance varies with GC content in a largely similar way for all constraint metrics, at least for low to intermediate GC content. Performance was defined to be the area under the precision-recall curve, auPRC, normalized by the fraction of examples that are positive, which varies among the lax truth sets. (C) We downsampled the lax truth set so that its size matched that of a stringent truth set, defined by essential genes (Methods), and used 1000 bootstraps of the truth sets to compute a mean and standard deviation of normalized auPRC for each constraint metric and each truth set (without stratifying by genomic features). (D) For each bootstrap sample, we computed the difference in normalized auPRC between bins of high and low values of GC content ((0.4, 0.7) and (0.2, 0.375), respectively), and reported the mean and standard deviation of those differences for each constraint metric relative to the stringent truth set. The results show that GC content is associated with changes in the performance of all constraint metrics relative to the stringent truth set. Code to reproduce C and D can be found here: https://github.com/quinlan-lab/constraint-tools/blob/main/papers/neutral_models_are_biased/11.compare-lax-with-stri ngent-truth-set.ipynb. In all panels, GC content was assessed in 1kb windows. Similar trends in performance are observed for the other constraint metrics, with the notable exception of CDTS and Depletion Rank, which increase in performance in GC-rich and/or gBGC-rich regions (Figure 4B, and Supporting Figure 5C and D, respectively). Interestingly, these data show that, of all the models, the Chen model has the best average performance, i.e., when evaluated at typical values of the features (vertical black lines in Figure 2A-C), even though it has the greatest bias (see in Table 2). Trading significant modeling bias on rare 12 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint windows for high model performance on common windows has been observed elsewhere (Kathail et al. 2024a). The number of windows in the truth set that we used to measure performance is large, enabling us to detect performance anomalies that occur deep in the tails of the feature distributions. But that truth set is also lax, as not all windows that overlap a Genehancer enhancer are expected to be highly constrained. For example, 18.4% of the noncoding genome is covered by Genehancer enhancers (Martini et al. 2023), but others have estimated that only 4.51% of the noncoding genome is under human-specific selection (Huber et al. 2020). We therefore identified noncoding windows that we conjectured would be under strong negative selection because they regulate the expression of essential genes (Methods). We combined them with an equal number of noncoding windows not overlapping any enhancer to form a small but stringent truth set. We used it to re-assess the performance of the various constraint metrics at predicting constraint. We observed that overall performance increases on moving from the lax to the stringent truth set, particularly for the Gnocchi and metrics (Figure 4C and Supporting λ𝑠 Figure 5E, F), providing a post-hoc validation of stringency. Even though the stringent truth set is small–making estimates of performance in the tails of the feature distributions noisy–we were nevertheless able to detect associations between performance and features (Figure 4D and Supporting Figure 5G,H), as we observed for the lax truth set (Figure 4B and Supporting Figure 5C,D).

Sequence conservation leverages the power of deep evolutionary time to provide strong signals of purifying selection, but this approach is limited to the discovery of selection common to a group of species. This forces us to look at human-specific diversity-based methods in order to find human-specific selection. The trade-off is that inference of selection via methods that rely on segregating sites is confounded in a number of ways. (1) Though rare, the fixation of an advantageous allele can drag linked, neutral DNA with it, sweeping away existing variation in that genomic region (Smith and Haigh 1974; Stephan 2019), and naively leading to false-positive calls of purifying selection. (2) During a population bottleneck, rare alleles tend to be lost quickly, before heterozygosity has been eroded, resulting in a transient imbalance between heterozygosity and allele number (Luikart et al. 1998; Amos and Hoffman 2010) that one could falsely interpret as the active removal of deleterious alleles in that interval by purifying selection, if demography is not accounted for. (3) The demographic history of a population can reduce the ability of natural selection to remove weakly deleterious variants (Lohmueller 2014), and so diversity-based methods, which rely on the removal of deleterious variants, will not discover constrained regions harboring those variants. Here, we identify additional factors that confound our interpretation of diversity-based methods for genetic-constraint calling. We show that current models of neutral genetic diversity are biased, compromising their ability to uncover constraint in the human noncoding genome in rare but important contexts, e.g., high GC content (enriched for promoters) or strong background selection (enriched for proximity to 13 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint genes), depending upon the model. This mirrors similar findings in the coding genome for a variety of organisms (Yap et al. 2010; Venkat et al. 2018; Ratnakumar et al. 2010). Similar to the difficulty in capturing edge cases when training machine-learning models to drive cars or answer questions in natural language, models of genetic diversity tend to capture the average behavior seen in a putatively neutral training set, at the expense of incurring large prediction errors on classes of genomic intervals that are rare in that training set. Capturing only average behavior is a blessing when rare events are true anomalies (e.g., constrained windows that are erroneously part of the training set), but a curse when they are not (e.g., truly neutral windows with high GC content). We find that the model with the highest performance on both the lax and stringent truth set—the Chen model—is also the model with the largest bias, suggesting that bias is a natural tradeoff of improved performance. Interestingly, the direct effect of GC content on bias is larger than that of gBGC and BGS, even though the former was explicitly included in the model whereas the latter were not. Since performance of the Chen model declines with GC content (Figure 4A), and since GC content is positively correlated with Gnocchi (Supporting Figure 6A,B), we expect that Gnocchi is most unreliable precisely where it predicts the most constraint. Moreover, given that GC content is auto correlated across the human genome, forming “isochores” (Supporting Figure 6B), we expect these false calls of constraint under the Chen model to cluster in genomic space.The relatively high performance of the Chen model (marginalized over all genomic feature values) may result from how comprehensive its set of features is, but it is likely that the way these features, particularly GC content, are included in the model is still not flexible enough to comprehensively capture neutral genetic diversity. Deep-learning models, such as those used to predict the density of noncoding SNVs (Fang et al. 2022), or their effects on fitness (Benegas et al. 2024; Trevino et al. 2021), do not suffer from such limitations on model complexity. Nevertheless, a deep-learning model mapping genomic features onto SNV density may still fall foul to the type of biases we illustrated here, as recently illustrated for other genomic prediction tasks (Kathail et al. 2024b; Bajwa et al. 2024). To the extent that this failure mode is due to insufficient drive to improve the fit in sparsely populated parts of feature space, a possible solution is to up-weight training examples there, e.g., by weighting a genomic interval in inverse proportion to the number of intervals similar to it (Morcos et al. 2011). In the meantime, we recommend annotating each genomic interval not only with its constraint score, but also a measure of how biased the score is expected to be in that interval, as a way to gauge confidence in the score. This can be done by grouping genomic intervals by the similarity of their feature vectors (with elements comprising GC content, etc), computing the distribution of constraint scores in each group, and assigning to each window the mean of the constraint-score distribution for the corresponding group. 14 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint

Provenance of constraint scores Gnocchi scores for a set of nonoverlapping 1000bp windows (“Chen windows”) were obtained from here (Supplementary Data #2). Windows on the X and Y chromosomes were omitted. Drs Nurdan Kuru and Adam Siepel ran ExtRaINSIGHT (Dukler et al. 2022) on the remaining Chen windows in order to compute , excluding Chen windows that did not pass the ExtRaINSIGHT λ𝑠 filters (private communication). A set of overlapping 500bp windows (“Halldorsson windows”) with Depletion Rank scores (Halldorsson et al. 2022) was obtained from here. A set of overlapping 550bp windows with CDTS scores (“CDTS windows”) (di Iulio et al. 2018) was obtained by private communication with Dr Amalio Telenti (Vir Biotechnology), and their coordinates were lifted over from hg19 to hg38. Assignment of genomic feature values to genomic windows We computed GC content in 1kb intervals centered on the Chen, Halldorsson and CDTS windows using “bedtools nuc”, e.g., as described here and here. We obtained BGS values computed by (Murphy et al. 2023) from here, which we converted to bed format and lifted over from hg19 to hg38 (see here). We obtained the gBGC coefficients on a set of genomic intervals for the European population computed by (Glémin et al. 2015) from here, and lifted over the genomic coordinates from hg18 to hg38, as described here. We assigned BGS values and gBGC coefficients (B_M1star.EUR) to the Chen, Halldorsson and CDTS windows by intersecting them with the bed files of BGS and gBGC, e.g., as shown here and here. Construction of the window sets to assess model bias, and the lax truth sets to assess constraint-prediction performance We filtered out untrustworthy windows, defined as windows that overlap any of the following intervals: gaps in the hg38 genome assembly, Encode “exclude regions” (Amemiya et al. 2019), and regions with insufficient read coverage in Gnomad version 3 (Chen et al. 2024). Exons were computed as shown here and merged using “bedtools merge”. Noncoding windows were defined to be Chen, Halldorsson and CDTS windows that don’t significantly overlap merged exons. Of the noncoding windows, those that don’t significantly overlap Genehancer enhancers (Fishilevich et al. 2017) were labeled “neutral”, and were used to assess model bias. Noncoding windows that do significantly overlap Genehancer enhancers were labeled “constrained”. Together the “neutral” and “constrained” windows comprised the lax truth set used to assess the performance of each of the four constraint metrics. Modeling the mapping from genomic features to constraint scores The code to compute the linear regressions of the constraint scores shown in the heat maps, including the regression coefficients and the proportion of variance explained, reported in the tables, can be found here. 15 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Computation of window residuals under the Chen model Each window has a residual, , equal to the difference between the expected SNV count under the Chen model of neutrality, , and the observed SNV count, . We computed from Gnomad Version 3, using the bash command “constraint-tools predict-germline-model-Nonly --windows ${chen_windows} REQUIRED_ARGUMENTS”, where REQUIRED_ARGUMENTS are specified here. We then used , in conjunction with Gnocchi, , to compute using the following algorithm: a = 1 b = -(2*S + G**2) c = S**2 s = sqrt(b**2 - 4*a*c) sign = 1 if G > 0 else -1 = (-b + sign*s)/(2*a) Finally: . Computation of the PR curve for the random classifier Given a truth set of size , let represent the proportion of windows that are constrained. Let represent the probability that a random classifier predicts that a given window is constrained. Then and where TP (average number of True Positives) is the average number of windows that are labeled “constrained” in the truth set and are predicted to be constrained by the classifier. Similar definitions exist for FP (average number of False Positives) and FN (average number of False Negatives). The formulae for Precision and Recall describe a horizontal line on the Precision-Recall curve (e.g., the dashed line in Figure 4A, traced out as the parameter is varied between 0 and 1) and represents a baseline level of performance that any meaningful constraint predictor ought to exceed. Computation of SNR The signal-to-noise ratio (SNR) is defined as , 16 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint where is the random residual of a constrained window, is the random residual of a neutral window, and represents the expectation operator. The numerator is the amplitude of the meaningful signal we are trying to recover (a dip in SNV count due to purifying selection); the denominator is the amplitude of the unwanted background noise we must contend with (unavoidable Poisson fluctuations in SNV count). Squaring the relation , where is the reduction in SNV count induced by purifying selection, we arrive at: . Applying the expectation operator, , we get: , where we used the facts that is constant (if is feature-dependent, then specialize this analysis to a particular bin of that feature), , and . Substituting these expressions for and into the definition of SNR, we get . When , SNR is approximately , indicating that the amplitudes of the signal and noise are approximately the same, making it impossible to recover the signal. This is reflected in a flat Precision-Recall (PR) curve (e.g., the dashed line in Figure 4A), which represents a baseline level of performance that any meaningful constraint predictor ought to exceed. We used the SNR formula to compute the numerical values of SNR reported in Figure 3 (and in Supporting Figure 4, where is feature-dependent). Why PR curves for all models overlap when there is little variability in the mock genomic feature (Figure 3G-I) Consider a set of simulated windows each with approximately the same value of the mock genomic feature. Under a given model, the predicted SNV count will be approximately the same for all these windows, because these models depend upon only that genomic feature. Now subset the windows into those that are constrained, and those that are not (neutral). The residual distributions for the two classes of windows (constrained and neutral) are just the observed SNV count distributions for the two classes, offset by a constant–the common predicted SNV count. Therefore the overlap of the residual distributions of the two classes is approximately the same for all models, and so all models perform approximately equally. Construction of the stringent truth set to assess performance of constraint scores Chen et al (Chen et al. 2024) obtained a set of highly constrained genes (ClinGen haploinsufficient genes, MGI essential genes, OMIM autosomal dominant genes, LOEUF 17 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint first-decile genes), linked enhancers to each gene using correlated patterns of histone modifications and gene expression, and assigned the most constrained enhancer to each gene, yielding a set of enhancers that we downloaded from here (Supplementary Data #6). We assigned genomic feature values, and constraint scores, to each of these “essential” enhancers similarly to how we did the assignments for the Chen, Halldorsson and CDTS windows. Finally, we constructed the stringent truth set by combining these “constrained” windows with an equal number of “neutral” windows sampled from the Chen lax truth set. See here. Funding This work was supported by a National Institutes of Health grant to ARQ (R01HG012252).

We graciously thank Nurdan Kuru and Adam Siepel for sharing their ExtRaINSIGHT data. We are grateful to Andrew Kern for expert feedback that significantly improved on an early draft of the manuscript.

Adrion, Jeffrey R., Christopher B. Cole, Noah Dukler, et al. 2020. “A Community-Maintained Standard Library of Population Genetic Models.” eLife 9 (e54967). https://doi.org/10.7554/eLife.54967. Amemiya, Haley M., Anshul Kundaje, and Alan P. Boyle. 2019. “The ENCODE Blacklist: Identification of Problematic Regions of the Genome.” Scientific Reports 9 (1): 9354. Amos, W., and J. I. Hoffman. 2010. “Evidence That Two Main Bottleneck Events Shaped Modern Human Genetic Diversity.” Proceedings. Biological Sciences 277 (1678): 131–137. Arndt, Peter F., Terence Hwa, and Dmitri A. Petrov. 2005. “Substantial Regional Variation in Substitution Rates in the Human Genome: Importance of GC Content, Gene Density, and Telomere-Specific Effects.” Journal of Molecular Evolution 60 (6): 748–763. Bajwa, Ayesha, Ruchir Rastogi, Pooja Kathail, Richard W. Shuai, and Nilah Ioannidis. 2024. “Characterizing Uncertainty in Predictions of Genomic Sequence-to-Activity Models.” Machine Learning in Computational Biology, March 15, 279–297. Benegas, Gonzalo, Carlos Albors, Alan J. Aw, Chengzhong Ye, and Yun S. Song. 2024. “GPN-MSA: An Alignment-Based DNA Language Model for Genome-Wide Variant Effect Prediction.” bioRxiv.org: The Preprint Server for Biology, ahead of print, April 6. https://doi.org/10.1101/2023.10.10.561776. Bishop, Christopher M. 2006. Pattern Recognition and Machine Learning. 1 (August): 740. Boffelli, Dario, Jon McAuliffe, Dmitriy Ovcharenko, et al. 2003. “Phylogenetic Shadowing of Primate Sequences to Find Functional Regions of the Human Genome.” Science (New York, N.Y.) 299 (5611): 1391–1394. 18 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Bulmer, M. 1986. “Neighboring Base Effects on Substitution Rates in Pseudogenes.” Molecular Biology and Evolution 3 (4): 322–329. Chandra, Sheel, and Ziyue Gao. 2025. “Sequence Context and Methylation Interact to Shape Germline Mutation Rate Variation at CpG Sites.” In bioRxiv. November 13. https://doi.org/10.1101/2025.11.13.688199. Charlesworth, Brian. 2013. “Background Selection 20 Years on: The Wilhelmine E. Key 2012 Invitational Lecture.” The Journal of Heredity 104 (2): 161–171. Chen, Siwei, Laurent C. Francioli, Julia K. Goodrich, et al. 2024. “A Genomic Mutational Constraint Map Using Variation in 76,156 Human Genomes.” Nature 625 (7993): 92–100. Cooper, Gregory M., Eric A. Stone, George Asimenos, et al. 2005. “Distribution and Intensity of Constraint in Mammalian Genomic Sequence.” Genome Research 15 (7): 901–913. Davydov, Eugene V., David L. Goode, Marina Sirota, Gregory M. Cooper, Arend Sidow, and Serafim Batzoglou. 2010. “Identifying a High Fraction of the Human Genome to Be under Selective Constraint Using GERP++.” PLoS Computational Biology 6 (12): e1001025. Dukler, Noah, Mehreen R. Mughal, Ritika Ramani, Yi-Fei Huang, and Adam Siepel. 2022. “Extreme Purifying Selection against Point Mutations in the Human Genome.” Nature Communications 13 (1): 4312. Duret, Laurent, and Nicolas Galtier. 2009. “Biased Gene Conversion and the Evolution of Mammalian Genomic Landscapes.” Annual Review of Genomics and Human Genetics 10: 285–311. Elango, Navin, Seong-Ho Kim, Eric Vigoda, and Soojin V. Yi. 2008. “Mutations of Different Molecular Origins Exhibit Contrasting Patterns of Regional Substitution Rate Variation.” PLoS Computational Biology 4 (2): e1000015. Fang, Yiyuan, Shuyi Deng, and Cai Li. 2022. “A Generalizable Deep Learning Framework for Inferring Fine-Scale Germline Mutation Rate Maps.” Nature Machine Intelligence 4 (12): 1209–1223. Fishilevich, Simon, Ron Nudel, Noa Rappaport, et al. 2017. “GeneHancer: Genome-Wide Integration of Enhancers and Target Genes in GeneCards.” Database: The Journal of Biological Databases and Curation 2017 (January): bax028. Friedman, Jerome H. 2001. “Greedy Function Approximation: A Gradient Boosting Machine.” Annals of Statistics 29 (5): 1189–1232. Fryxell, Karl J., and Won-Jong Moon. 2005. “CpG Mutation Rates in the Human Genome Are Highly Dependent on Local GC Content.” Molecular Biology and Evolution 22 (3): 650–658. Glémin, Sylvain, Peter F. Arndt, Philipp W. Messer, Dmitri Petrov, Nicolas Galtier, and Laurent Duret. 2015. “Quantification of GC-Biased Gene Conversion in the Human Genome.” Genome Research 25 (8): 1215–1228. Gower, Graham, Nathaniel S. Pope, Murillo F. Rodrigues, et al. 2025. “Accessible, Realistic Genome Simulation with Selection Using Stdpopsim.” In bioRxivorg. August 12. https://doi.org/10.1101/2025.03.23.644823. 19 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Halldorsson, Bjarni V., Hannes P. Eggertsson, Kristjan H. S. Moore, et al. 2022. “The Sequences of 150,119 Genomes in the UK Biobank.” Nature 607 (7920): 732–740. Huber, Christian D., Bernard Y. Kim, and Kirk E. Lohmueller. 2020. “Population Genetic Models of GERP Scores Suggest Pervasive Turnover of Constrained Sites across Mammalian Evolution.” PLoS Genetics 16 (5): e1008827. Iulio, Julia di, Istvan Bartha, Emily H. M. Wong, et al. 2018. “The Human Noncoding Genome Defined by Genetic Diversity.” Nature Genetics 50 (3): 333–337. Kathail, Pooja, Richard W. Shuai, Ryan Chung, Chun Jimmie Ye, Gabriel B. Loeb, and Nilah M. Ioannidis. 2024a. “Current Genomic Deep Learning Models Display Decreased Performance in Cell Type Specific Accessible Regions.” bioRxiv.org: The Preprint Server for Biology, July 10, 2024.07.05.602265. Kathail, Pooja, Richard W. Shuai, Ryan Chung, Chun Jimmie Ye, Gabriel B. Loeb, and Nilah M. Ioannidis. 2024b. “Current Genomic Deep Learning Models Display Decreased Performance in Cell Type-Specific Accessible Regions.” Genome Biology 25 (1): 202. Lauterbur, M. Elise, Maria Izabel A. Cavassim, Ariella L. Gladstein, et al. 2023. “Expanding the Stdpopsim Species Catalog, and Lessons Learned for Realistic Genome Simulations.” May 23. https://doi.org/10.7554/elife.84874. Liu, Li, Maxwell D. Sanderford, Ravi Patel, Pramod Chandrashekar, Greg Gibson, and Sudhir Kumar. 2019. “Biological Relevance of Computationally Predicted Pathogenicity of Noncoding Variants.” Nature Communications 10 (1): 330. Lohmueller, Kirk E. 2014. “The Distribution of Deleterious Genetic Variation in Human Populations.” Current Opinion in Genetics & Development 29 (December): 139–146. Luikart, G., F. W. Allendorf, J. M. Cornuet, and W. B. Sherwin. 1998. “Distortion of Allele Frequency Distributions Provides a Test for Recent Population Bottlenecks.” The Journal of Heredity 89 (3): 238–247. Martini, Lorenzo, Roberta Bardini, Alessandro Savino, and Stefano Di Carlo. 2023. “GRAIGH: Gene Regulation Accessibility Integrating GeneHancer Database.” 2023 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), December 5, 343–348. McVicker, Graham, David Gordon, Colleen Davis, and Phil Green. 2009. “Widespread Genomic Signatures of Natural Selection in Hominid Evolution.” PLoS Genetics 5 (5): e1000471. Meader, Stephen, Chris P. Ponting, and Gerton Lunter. 2010. “Massive Turnover of Functional Sequence in Human and Other Mammalian Genomes.” Genome Research 20 (10): 1335–1343. Morcos, Faruck, Andrea Pagnani, Bryan Lunt, et al. 2011. “Direct-Coupling Analysis of Residue Coevolution Captures Native Contacts across Many Protein Families.” Proceedings of the National Academy of Sciences of the United States of America 108 (49): E1293–301. Mugal, Carina F., and Hans Ellegren. 2011. “Substitution Rate Variation at Human CpG Sites Correlates with Non-CpG Divergence, Methylation Level and GC Content.” Genome Biology 12 (6): R58. 20 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Murphy, David A., Eyal Elyashiv, Guy Amster, and Guy Sella. 2023. “Broad-Scale Variation in Human Genetic Diversity Levels Is Predicted by Purifying Selection on Coding and Non-Coding Elements.” eLife 12 (June). https://doi.org/10.7554/eLife.76065. Pollard, Katherine S., Melissa J. Hubisz, Kate R. Rosenbloom, and Adam Siepel. 2010. “Detection of Nonneutral Substitution Rates on Mammalian Phylogenies.” Genome Research 20 (1): 110–121. Ponting, Chris P., and Ross C. Hardison. 2011. “What Fraction of the Human Genome Is Functional?” Genome Research 21 (11): 1769–1776. Pouyet, Fanny, Simon Aeschbacher, Alexandre Thiéry, and Laurent Excoffier. 2018. “Background Selection and Biased Gene Conversion Affect More than 95% of the Human Genome and Bias Demographic Inferences.” eLife 7 (August). https://doi.org/10.7554/eLife.36317. Rands, Chris M., Stephen Meader, Chris P. Ponting, and Gerton Lunter. 2014. “8.2% of the Human Genome Is Constrained: Variation in Rates of Turnover across Functional Element Classes in the Human Lineage.” PLoS Genetics 10 (7): e1004525. Ratnakumar, Abhirami, Sylvain Mousset, Sylvain Glémin, et al. 2010. “Detecting Positive Selection within Genomes: The Problem of Biased Gene Conversion.” Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 365 (1552): 2571–2580. Schraiber, Joshua G., Jeffrey P. Spence, and Michael D. Edge. 2024. “Estimation of Demography and Mutation Rates from One Million Haploid Genomes.” In Genomics, No. Biorxiv;2024.09.18.613708v1. BioRxiv, September 22. https://www.biorxiv.org/content/10.1101/2024.09.18.613708v1.full.pdf. Schrider, Daniel R., and Andrew D. Kern. 2015. “Inferring Selective Constraint from Population Genomic Data Suggests Recent Regulatory Turnover in the Human Brain.” Genome Biology and Evolution 7 (12): 3511–3528. Smith, J. M., and J. Haigh. 1974. “The Hitch-Hiking Effect of a Favourable Gene.” Genetical Research 23 (1): 23–35. Smith, Nick G. C., Mikael Brandström, and Hans Ellegren. 2004. “Evidence for Turnover of Functional Noncoding DNA in Mammalian Genome Evolution.” Genomics 84 (5): 806–813. Stephan, Wolfgang. 2019. “Selective Sweeps.” Genetics 211 (1): 5–13. Sved, J., and A. Bird. 1990. “The Expected Equilibrium of the CpG Dinucleotide in Vertebrate Genomes under a Mutation Model.” Proceedings of the National Academy of Sciences of the United States of America 87 (12): 4692–4696. Telis, Natalie, Robin Aguilar, and Kelley Harris. 2020. “Selection against Archaic Hominin Genetic Variation in Regulatory Regions.” Nature Ecology & Evolution 4 (11): 1558–1566. Trevino, Alexandro E., Fabian Müller, Jimena Andersen, et al. 2021. “Chromatin and Gene-Regulatory Dynamics of the Developing Human Cerebral Cortex at Single-Cell Resolution.” Cell 184 (19): 5053–5069.e23. 21 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint Venkat, Aarti, Matthew W. Hahn, and Joseph W. Thornton. 2018. “Multinucleotide Mutations Cause False Inferences of Lineage-Specific Positive Selection.” Nature Ecology & Evolution 2 (8): 1280–1288. Ward, Lucas D., and Manolis Kellis. 2012. “Evidence of Abundant Purifying Selection in Humans for Recently Acquired Regulatory Functions.” Science (New York, N.Y.) 337 (6102): 1675–1678. Yap, Von Bing, Helen Lindsay, Simon Easteal, and Gavin Huttley. 2010. “Estimates of the Effect of Natural Selection on Protein-Coding Content.” Molecular Biology and Evolution 27 (3): 726–734. 22 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 28, 2026. ; https://doi.org/10.64898/2026.01.28.701168doi: bioRxiv preprint