iiinova — where your digital self keeps evolving

Jeff Spence @jeffspence.github.io · 22/06/2026

Finally, a huge thank you to @jonj-udd.bsky.social for leading this project, @jkpritch.bsky.social for supervising, and our coauthors @nikhilmilind.dev, @lindakachuri.bsky.social, and @johnwitte.bsky.social. Thanks also to everyone who pointed us to references or provided feedback or data.

The acknowledgments paragraph from the paper:

We thank Josh Schraiber for helpful discussions, Molly Przeworski for helpful comments on the manuscript, and Shai Carmi, Graham Coop, and Hilary Martin for pointing us to related work. We additionally thank our anonymous reviewers for their comments that have helped to improve this work. Anastasia Stolyarova kindly provided processed data from another study (Stolyarova et al. 2025). Thanks to the Pritchard Lab for many valuable comments and discussions.

190

Jeff Spence @jeffspence.github.io · 22/06/2026

Nevertheless, we find that if a population is already at equilibrium and then experiences demographic changes, those changes don't affect the expected frequency too much. For example, the expected allele frequency after 176 generations of exponential growth is close to the initial frequency.

Plot with the expected allele frequency on the y-axis plotted against the strength of selection against homozygotes on the x-axis. There are four lines that are all quite close: one line assuming equilibrium for the initial size of 17,000, one line assuming equilibrium fro the final size of 508,863, one line for a population that grows exponentially from 17,000 to 508,863 over 176 generation, and one line assuming an infinite population size at equilibrium. The exponential growth line largely overlaps the line for the initial size. The line for the final size largely overlaps the line for the infinite population size.

110

Jeff Spence @jeffspence.github.io · 22/06/2026

For purely recessive selection, we show that equilibration takes a bit longer. For a recessive lethal gene with a typical mutation rate, equilibration takes hundreds to a few thousand generations.

Chart with expected allele frequency on the y-axis, and the number of generations following the onset of selection on the x-axis. Different lines show different strengths of homozygous selection and different initial allele frequencies. Expected allele frequencies generally begin converging in tens to hundreds of generations, and fully converging in hundreds to thousands of generations.

Chart with the KL divergence between the distribution of allele frequencies on the y-axis and the number of generations following the onset of selection on the x-axis. Again, different lines show different strengths of purely recessive selection and different initial allele frequencies. The KL divergence tends to start decreasing after tens to hundreds of generations.

130

Jeff Spence @jeffspence.github.io · 31/03/2026

We then forecasted how much genetic diversity will be lost across species given current estimates of habitat loss. In the short term, we predict very little will be lost, but in the long term, species might lose up to 20% of their genetic diversity on average. This is WITHOUT FURTHER HABITAT LOSS!

Forecasts showing how genetic diversity will change in the short and medium terms. In the short term 0% to 14% of genetic diversity will be lost, but in the mid-term up to 20% will be lost.

152

Jeff Spence @jeffspence.github.io · 31/03/2026

What is happening is that locally the species is actually becoming significantly _less_ genetically diverse, but the fragmentation disconnects parts of the species, which slowly diverge from each other. The species-wide diversity masks the plummeting local diversity. Again, this is BAD.

Plots showing that within population fragments genetic diversity plummets following habitat fragmentation.

132

Jeff Spence @jeffspence.github.io · 31/03/2026

In these fragmentation scenarios, we actually see genetic diversity _increase_ in the long term. This is super weird! How can destroying habitat result in increased genetic diversity?

Plots showing that genetic diversity increases following habitat fragmentation.

110

Jeff Spence @jeffspence.github.io · 31/03/2026

Everything up to now has assumed that habitat is lost along one edge of the species range (as a rough model for the effects of climate change). If instead we model the kinds of habitat fragmentation that result from land use changes, things get weird.

Schematics showing species losing habitat in random fragments resulting in a "patchy" landscape.

140

Jeff Spence @jeffspence.github.io · 31/03/2026

What we found surprised us... Species can lose A LOT of habitat without having much effect on genetic diversity. In some scenarios a species can lose 80% of its habitat while retaining 90% of its genetic diversity. But this is just IN THE SHORT TERM.

Two plots. Left plot shows genetic diversity on the y-axis against time on the x-axis. It takes thousands of generations following habitat loss for the full effects to be felt on genetic diversity. Right plot shows genetic diversity on the y-axis against habitat loss on the x-axis. In the long term habitat loss results in an almost one-to-one reduction in genetic diversity, whereas losses are much less extreme in the short term.

121

Jeff Spence @jeffspence.github.io · 31/03/2026

We started with a simple question: What happens to genetic diversity when part of a species' habitat suddenly become uninhabitable?

Schematic showing a species losing part of its habitat over time.

130

Jeff Spence @jeffspence.github.io · 13/12/2025

We also get a nice formula in the no drift limit (huge population sizes).

$An equation rendered from the following LaTeX: \lim_{N\to\infty} \mathbb{E}[f] = -\frac{\varphi}{2(1-\varphi)} + \sqrt{\left(\frac{\varphi}{2(1-\varphi)}\right)^2 + \frac{\mu}{(1-\varphi)s_{hom}}}$

120

Jeff Spence @jeffspence.github.io · 13/12/2025

In the limit as φ goes to zero, this reduces to a result of Nei that @jgschraiber.bsky.social graciously pointed out to us. www.pnas.org/doi/abs/10.1...

$an equation rendered from the following LaTeX: \lim_{\varphi \downarrow 0} \mathbb{E}[f] = \frac{\Gamma\left(2N\mu + 1/2\right)}{\Gamma\left(2N\mu\right)\sqrt{2Ns_{hom}}}$

120

Jeff Spence @jeffspence.github.io · 13/12/2025

We provide some approximations for the expected frequency under mutation-selection balance in this setup. The first is a convergent series that's accurate for arbitrary φ and strong s. It's admittedly a bit gnarly (aesthetically and numerically).

$A gnarly formula showing the expected value of the allele frequency. It is rendered from the following LaTeX: \mathbb{E}[f] = \frac{\sum_{k=0}^\infty \frac{\left(-4Ns_{hom}\varphi/(1+\varphi\right)^k \Gamma\left(2N\mu/(1+\varphi) + (k+1)/2\right)}{k! \left(2Ns_{hom}(1-\varphi)/(1+\varphi)\right)^{2N\mu/(1+\varphi) + (k+1)/2}}}{\sum_{\ell=0}^\infty \frac{\left(-4Ns_{hom}\varphi/(1+\varphi\right)^\ell \Gamma\left(2N\mu/(1+\varphi) + \ell/2\right)}{\ell! \left(2Ns_{hom}(1-\varphi)/(1+\varphi)\right)^{2N\mu/(1+\varphi) + \ell/2}}}$

140

Jeff Spence @jeffspence.github.io · 07/11/2025

Basically, when a variant gets lucky and drifts to high frequency, it is more likely to be a hit for every trait that it affects. This makes top hits APPEAR more pleiotropic, when in fact they’re actually LESS pleiotropic on average. 18/n

Three plots showing that variation in minor allele frequency can result in the top GWAS hits being hits for more traits on average. The leftmost plot shows that minor allele frequency decreases for less significant GWAS hits in both simulations and empirically. The middle plot shows that in simulations the most significant GWAS hits are, on average, the most trait specific. The right plot shows that in simulations and empirically the top GWAS hits are significant for more traits than less significant GWAS hits.

150

Jeff Spence @jeffspence.github.io · 07/11/2025

This randomness in variant frequencies also results in what seems a paradox: top GWAS hits are both more trait specific AND more likely to be hits for other traits! 17/n

A schematic showing the p-values for two traits. Specificity makes p-values more significant for one trait at the expense of the other trait. Increasing minor allele frequency makes a variant more significant for both traits.

140

Jeff Spence @jeffspence.github.io · 07/11/2025

GWAS prioritize individual variants, and sometimes variants just happen to get lucky and drift to high frequencies. This plays a surprisingly large role in prioritization in simulated GWAS. 16/n

Two plots. The left shows trajectories of identical mutations that all arose at the same time in the past, and how they evolve forward in time to the present day. Despite being identical mutations, their trajectories show a wide spread with a large amount of variance in their frequencies at the present day. The right shows the "realized heritability" of simulated SNPs and how it relates to the squared effect on the study trait. There is little concordance between realized heritability and squared effect on the study trait. Instead, there is a strong relationship between realized heritability and minor allele frequency, which, in the plot is denoted using color.

160

Jeff Spence @jeffspence.github.io · 07/11/2025

Specificity is not the end of the story. Burden tests and GWAS also prioritize genes based on things that have absolutely nothing to do with traits. Burden tests aggregate signal across variants. Long genes have more variants, and so tend to get prioritized higher. 15/n

Three scatter plots showing that long genes do not have larger effects on traits on average, but they do have smaller standard errors, and as a result they have stronger squared z scores.

141

Jeff Spence @jeffspence.github.io · 07/11/2025

In contrast, HHIP (hedgehog interacting protein) has tons of GWAS hits near it, and basically no burden signal. HHIP makes sense as a height hit, but it’s also been implicated in COPD, pituitary hormone deficiency, etc… It’s NOT height specific! 14/n

A zoomed-in version of the HHIP locus from a previous post. There are many genome-wide significant GWAS hits, some with p-values smaller than 10^-150. No genes in the locus show significant burden signal.

120

Jeff Spence @jeffspence.github.io · 07/11/2025

NPR2 is the second most significant gene in the burden tests, but is in only the 243rd most significant GWAS locus. Homozygous LoFs in NPR2 cause severe short stature, but don’t affect intelligence, facial features, etc… seems like NPR2 is height specific! 13/n

A zoomed in version of the NPR2 locus plot from the previous post. There are two unlinked GWAS hits that are just barely genome-wide significant, and NPR2 shows strong burden test signal.

140

Jeff Spence @jeffspence.github.io · 07/11/2025

To recap: Burden tests prioritize trait-SPECIFIC GENES. GWAS prioritize genes near trait-SPECIFIC VARIANTS. Looking at height gives a couple of really nice examples: NPR2 and HHIP. 12/n

A scatter plot showing p-values for GWAS loci on the x-axis and p-values for burden tests on the y-axis. There is some concordance, but a lot of scatter. Two points are circled and point to additional panels. One point is labeled NPR2 and points to plots showing a locus where there is very little GWAS signal (only two unlinked variants that are just barely genome-wide significant), but the NPR2 gene shows strong association in burden tests. Another point is labeled HHIP and points to plots showing a locus where there are many extremely significant GWAS hits, but essentially no burden signal.

141

Jeff Spence @jeffspence.github.io · 07/11/2025

Both of these ways of being specific contribute to variants being ranked highly in GWAS: 1. Coding variants in specifically-expressed genes are more highly ranked AND 2. Non-coding variants in tissue-specific ATAC peaks are more highly ranked. 11/n

Two LDSC plots showing that heritability enrichment increases with context specificity. The left plot shows that coding variants contribute more to heritability if they are in more specifically-expressed genes. The right plot shows that non-coding variants contribute more to heritability if they are in tissue-specific ATAC peaks.

140

Jeff Spence @jeffspence.github.io · 07/11/2025

What about GWAS? GWAS prioritize genes near trait-specific VARIANTS, whereas burden tests prioritize trait specific GENES. Variants can be specific because they act on trait-specific genes, or because they act on pleiotropic genes in a context-specific way. 10/n

A schematic showing the two ways that variants can be trait specific. They can either affect a pleiotropic gene in a context-specific way, or they can affect a trait-specific gene.

160

Jeff Spence @jeffspence.github.io · 07/11/2025

Using theory developed by @yuvalsim.bsky.social and @gs2747.bsky.social (journals.plos.org/plosbiology/...) we predicted that burden tests rank genes by SPECIFICITY! This is surprising! Burden tests DO NOT rank genes by IMPORTANCE! These predictions play out in the UKB. 8/n

Left side shows a cartoon of the derivation for why burden test rank genes in part by trait specificity. Selection decreases LoF frequencies, with frequency being inversely proportional to the selection coefficient. Selection is driven by importance across all traits, with the selection coefficient being proportional to the sum of squared effects on all traits. The right hand side shows a stratified QQ plot, where more specifically-expressed genes tend to have much stronger burden test signals.

171

Jeff Spence @jeffspence.github.io · 07/11/2025

Instead, one might want to rank genes by their trait SPECIFICITY -- how much do they affect the study trait relative to their effects across all traits. A trait-specific gene might be more "core" to trait biology. 7/n

A cartoon showing the mathematical definitions of "trait importance" (the squared effect on the trait under study) and "trait specificity" (the squared effect on the trait under study divided by the sum of squared effects across all traits).

150

Jeff Spence @jeffspence.github.io · 07/11/2025

We checked for ourselves, and found: 1. Most burden hits are near a GWAS hit (they converge!) BUT 2. The ranking of hits is surprisingly discordant. E.g., the second most significant burden hit for height is ranked 243rd in GWAS! 4/n

A bar chart showing burden hits for 12 traits. Each bar represents an individual gene that was found to be significant for a given trait using LoF burden tests, and they are stacked in order of significance. They are colored by the rank of the locus that contains them (if any) in GWAS. There are many top burden hits that are not ranked highly in GWAS, with some top burden hits not even being in the top 200 or 300 most significant GWAS loci.

131

Jeff Spence @jeffspence.github.io · 07/11/2025

GWAS and burden tests both regress trait values against genetic variation. In line with this conceptual similarity, previous work (link.springer.com/article/10.1... , www.biorxiv.org/content/10.1...) suggested these tests “converge” on the same genes. 2/n

Schematic showing how GWAS and burden tests work. Both regress phenotypes against genetic variation. In GWAS "genetic variation" is genotypes at individual variants. In burden tests individuals are grouped together depending on whether or not they have any LoF within a given gene.

130

Jeff Spence @jeffspence.github.io · 30/05/2025

Overall, many of our models did a pretty job of ranking which guides would drive expression the most, but almost all of that performance came across genes. When trying to predict which guide would have the largest effect on the expression of a particular gene, our results were more mixed.

A scatter plot of the predicted fold-change in expression of different genes for different guide RNAs compared to the actual fold-change measured by qPCR. There is generally a good concordance between the two (Spearman correlation of 0.7956), but the points are colored by which gene they came from, and it is clear that much of the predictive power comes from cross-gene accuracy, whereas accuracy for predicting the effects of different guides acting on the same gene is lower.

100

Jeff Spence @jeffspence.github.io · 30/05/2025

And we considered a lot of different models of what dCas9-p300 actually does to chromatin tracks (How far away does p300 acetylate? How strongly would it increase signal in a ChIP-seq track? Does this interact with nucleosome occupancy? Do nucleosomes interfere with guide binding?)

Schematic of models of how we model the effects of dCas9-p300 on local chromatin. Panel A shows a dCas9-p300 preferentially binding to a region with lower nucleosome occupancy (as measured by MNase-seq). Panel B shows dCas9-p300 depositing H3K27ac marks with some intensity lambda, and spread over a distance sigma. Panel C shows the resulting perturbed H3K27ac track.

100

Jeff Spence @jeffspence.github.io · 30/05/2025

On the experimental side, we can use dCas9 fused to a chromatin modifier to locally alter chromatin structure. We can then read out how those modifications affect gene expression. This lets us directly test the “causal understanding” of our deep learning models!

Cartoon showing an overview of the study.
The top row shows a flow chart of how epigenetic data is used to train neural networks to predict either endogenous expression or predict the effects of how epigenome editing affects expression.
The middle row shows that the neural networks take histone modification epigenetic tracks as input and then predict expression levels.
The bottom row shows how perturbed epigenetic tracks get used to make predictions that can then be compared to experimental results obtained using dCas9-p300.

120

Jeff Spence @jeffspence.github.io · 28/03/2025

E.g., for height we also don't find a lot of loci that have significant burden hits but lack a GWAS hit. But we do find the opposite quite a bit. We were also struck by how discordant the rankings are, and arguably when you have this many significant loci, some kind of ranking is necessary. 5/6

Scatter plot showing minimum p-values in LD blocks for GWAS and burden tests. Very few blocks are genome-wide significant in the burden test but not genome-wide significant in the GWAS, whereas many blocks are genome-wide significant in GWAS but not in burden tests.

141

Jeff Spence @jeffspence.github.io · 17/12/2024

The plots. In all plots the x-axis is "p-value rank".

Left: The y-axis is "Mean MAF relative to overall mean MAF" and the plot shows both simulations and empirical results. The simulations and empirical results show good agreement, and both show that the most significant GWAS hits are are higher than average MAF.

Center: The y-axis is "Mean specificity" and the plot only shows simulations. The most significant GWAS hits are the most specific.

Right: The y-axis is "Mean number of significant traits per hit". The plot shows both simulations and empirical results and they are in good agreement. Both show that the top GWAS hits are hits for more than 2 traits on average, while weaker GWAS hits are only hits for between 1.5 and 1.75 traits on average.

132

Jeff Spence @jeffspence.github.io · 17/12/2024

This randomness in variant frequencies also results in what seems a paradox: top GWAS hits are both more trait specific AND more likely to be hits for other traits! 15/n

A schematic showing two traits. Increasing trait specificity for the first trait increases power for that trait but decreases power for the other trait. Increasing MAF increases power for both traits.

120

Jeff Spence @jeffspence.github.io · 17/12/2024

GWAS prioritize individual variants, and sometimes variants just happen to get lucky and drift to high frequencies. This plays a surprisingly large role in prioritization in simulated GWAS 14/n

A scatter plot titled "In GWAS, sufficiently strong variants are mainly prioritized by frequency, no effect size". The y-axis is "realized heritability", a measure of how highly a variant would be prioritized, and the x-axis is squared effect size on the study trait. Individual points are simulated SNPs. The most highly ranked SNPs do no necessarily have the highest squared effect sizes, but do seem to have the highest MAF.

120

Jeff Spence @jeffspence.github.io · 17/12/2024

Three plots.

Left: a scatter plot of an unbiased estimate of mean squared effect size across traits within a bin (y-axis) plotted against the mean expected LoFs in the bin (x-axis). An arrow indicates that the x-axis also corresponds to genes with longer CDS. There is essentially no trend between the two.

Center: a scatter plot with the same x-axis, but the y-axis is now the mean standard error of LoF burden estimates of the effect size. There is a strong negative trend with longer genes having smaller standard errors.

Right: a scatter plot with the same x-axis, but the y-axis is now the mean squared z-score across traits. There is a strong positive relationship with longer genes having more burden signal.

130

Jeff Spence @jeffspence.github.io · 17/12/2024

To recap so far: burden tests prioritize trait-specific genes; GWAS prioritize genes near trait-specific variants. Height GWAS and burden tests give a couple of really nice examples: NPR2 and HHIP. 10/n

Three plots:

Top: The same scatter plot as before showing how burden tests and GWAS prioritize different loci for height. Two genes, HHIP and NPR 2 are highlighted.

Bottom left: a closer look at the NPR2 locus. There are 2 weak unlinked GWAS hits in the locus, but tons of burden signal for the NPR2 gene, and not much burden signal at other genes in the locus.

Bottom right: a closer look at the HHIP locus. There are tons of very strong unlinked GWAS hits, but none of the genes in the locus show much if any burden signal.

230

Jeff Spence @jeffspence.github.io · 17/12/2024

We find that both of these ways of being specific contribute to variants being ranked highly in GWAS: 1) Coding variants in specifically-expressed genes are more highly ranked 2) Non-coding variants in tissue-specific ATAC peaks are more highly ranked 9/n

Two plots.

Left side: a plot showing LDSC's tau plotted against expression specificity bins for 9 trait-tissue pairs. Across all traits, and on average, tau increases with increasing tissue specificity.

Right side: a similar plot showing LDSC's tau plotted against ATAC-seq peak specificity for non-coding variants. Across all traits, and on average, tau increases with increasing ATAC peak specificity.

130

Jeff Spence @jeffspence.github.io · 17/12/2024

What about GWAS? GWAS prioritize genes near trait-specific VARIANTS. This is profoundly different from prioritizing trait-specific GENES. Variants can be specific because they act on trait-specific genes, or because they act on pleiotropic genes in a context-specific way. 8/n

A schematic titled "Variants can be trait specific in two ways". The left hand side shows contours of variant specificity in a plot with the x-axis being the trait specificity of the gene through which the variant acts and the y-axis being how specific the variant is relative to that gene. Three points are labeled, and presented with separate cartoons. Points 1 and 2 show variants acting through a highly pleiotropic gene. Variant 2 is coding and hence also pleiotropic and not trait specific. Variant 1 is non-coding and only affects a specific context resulting in it being quite trait specific. Variant 3 is a coding variant but in a gene that is itself trait specific, making the variant also trait specific.

140

Jeff Spence @jeffspence.github.io · 17/12/2024

We developed theory to understand what determines burden test gene ranking: surprisingly, these rank by SPECIFICITY! We confirmed this in practice in UKB. This is surprising! Burden tests don’t rank genes by their importance! 6/n

Left side: A cartoon showing the derivation of a theoretical result. Burden signal depends on squared effect size on the focal trait and the frequency of LoFs. Selection decreases LoF frequency, making frequency inversely related to the selection coefficient. But selection is driven by a variant's effects on traits, so the selection coefficient is proportional to the sum of squared effects across all traits. This results in burden signal being proportional to the squared effect on the focal trait relative to the sum of squared effects across all traits. We call this "trait specificity".

Right side: a quantile-quantile plot of burden test results stratified by expression specificity bin. The signal is much stronger in the most specifically-expressed bins. Some genes are labeled (e.g., ACAN for height, HBB for MCV, and GPT for ALT).

150

Jeff Spence @jeffspence.github.io · 17/12/2024

Maybe instead, one might want to rank genes by their trait specificity – how important they are for the study trait relative to how important they are across all traits. A trait-specific gene might better capture how “core” it is to trait biology. 5/n

Left side: The title says "Which gene should be prioritized?" and below is a cartoon showing two genes. One gene only affects the study trait. The other gene has a much stronger effect on the study trait but also affects many other genes.

Right side: an equation with the definitions of trait importance (the squared effect size on the study trait) and trait specificity (the ratio of the squared effect size on the study trait to the sum of squared effect sizes across all fitness-relevant traits).

280

Jeff Spence @jeffspence.github.io · 17/12/2024

GWAS and burden tests are conceptually similar, so we might expect them to prioritize genes similarly. We find that their prioritizations of genes are surprisingly discordant (even when trying to get them to match as much as possible!) 2/n

Two plots. The first plot is a bar chart showing top genes as ranked by burden tests and whether or not they are in one of the top GWAS loci. There are several bars for a number of traits including height, MCH, LDL, and FVC. Generally less than half of top burden genes are also in top GWAS loci. The second plot is a scatter plot comparing minimum p-values for a height GWAS in different GWAS loci against minimum p-values from LoF burden tests for those loci. The points are somewhat correlated, but there are many points that are genome-wide significant in the GWAS but not in the burden test. There are also some strongly significant loci by the burden test that are much more lowly ranked by GWAS.

161

Jeff Spence @jeffspence.github.io · 17/12/2024

But this means that what matters for association study ranking is how specific a variant’s effects are to the trait you’re studying! signal ~ μ γ² / Σₜ γₜ² The most important variants don't end up with the largest z-scores! 5/6

A scatter plot of mean z-score squared (y-axis) plotted against an unbiased estimate of mean squared effect size (x-axis). The trend is essentially totally flat for genes with large effects.

131

Jeff Spence @jeffspence.github.io · 17/12/2024

But selection only cares about a variant because of how it impacts fitness-relevant traits. Under certain assumptions (see the paper!), we have that the strength of selection is the sum of effects across all fitness-relevant traits: s = Σₜ γₜ² www.biorxiv.org/content/10.1... 4/6

A scatter plot showing the mean of an unbiased estimate of squared effect sizes across traits for genes in bins (y-axis) plotted against the mean s_het within the bin (x-axis). The plot shows a strong positive trend between mean squared effect size and strength of selection.

130

Jeff Spence @jeffspence.github.io · 17/12/2024

If a variant experiences strong selection, with selection coefficient s then p(1-p) ~ μ/s , where μ is the mutation rate (thanks JBS Haldane!!!). We see this in the UKB! This means: signal ~ μγ²/s 3/6

Scatter plot showing mean LoF frequency in bins (y-axis) plotted against mean s_het in bins (x-axis). The plot shows a strong negative relationship between frequency and strength of selection.

140