flowchart TD A[TE genotype matrix + flanking SNPs<br/>per-window image] --> B[Two-branch CNN<br/>state channel + continuous channel] B --> C[Per-window embeddings] C --> D[attention pooling/deepsets<br/>variable number of TEs K] D --> E[Neural posterior estimator<br/>normalising flow, sbi] E --> F[DFE parameters<br/>gamma shape, mean s,<br/>burst timing/duration/amplitude,<br/>insertion rate]
Inferring the fitness effects and dynamics of transposable elements with simulation-based deep learning
An informal progress report
Overview
A long-standing goal of my work has been to obtain more precise estimates of the fitness effects and dynamics of transposable elements (TEs). TEs are known to be eco-responsive, to be at least mildly deleterious, to undergo bursts of transposition and to display insertion preferences. These properties make them substantially harder to model than SNPs. A promising avenue of research probably consists in contrasting patterns of TE age, frequency, and genome-wide distribution to obtain quantitative estimates of the parameters above.
This informal post summarises a suite of tools I have been developing towards that goal, and the directions I am currently pursuing. It is perhaps a bit naive, and I still need it to be reviewed by experts in the field of deep learning and population genetics, but the first results are encouraging enough from my perspective that I think they may be shared with the community. I hope this can stimulate more, better work on this, and maybe illustrate the challenges and excitement of taking advantage of some of the most powerful methods in population genomics for what remains a “niche” question, the population dynamics of TEs.
Potential of Simulation-Based Inference for the population dynamics of TEs.
A workhorse of population genomics is Approximate Bayesian Computation (ABC): by comparing observed summary statistics with simulations, one retains the simulations closest to the data and recovers posterior distributions for the parameters of interest. This works well for SNP data. TEs are more delicate, for three reasons.
First, TEs are far less abundant than SNPs along the genome. If one wishes to simulate TEs in their genomic context as realistically as possible, it may require simulating megabases of host’s sequence to obtain tens of polymorphic TEs for a given TE lineage. Second, the insertion rate can vary through time, for instance during a burst of transposition, which produces cohorts of similarly aged elements. A recent burst can generate an excess of singletons in the allele frequency spectrum that has nothing to do with selection, yet would ordinarily be read as such (Bourgeois and Boissinot 2019). Third, simulating selection against TEs requires forward-in-time simulators such as SLiM, which are slow and computationally costly (cf infra).
So far I had only compared simulations with observations in a qualitative or semi-quantitative way (Bourgeois et al. 2020; Horvath et al. 2024). That past work suggests that the abundance of TEs along the genome, together with the contrast between their frequency and their age, should carry enough signal to disentangle transposition activity, preferential insertion, and selection. An age-of-allele test of neutrality has also been proposed a few years ago already (Horvath et al. 2022), and implemented in a Bayesian framework that showed promise (Blumenstiel et al. 2014). Using ABC (or any simulation-based inference methodology) to compare TE simulations with observations should therefore be doable to jointly infer selection, insertion bias, and time-varying transposition dynamics.
However, ABC requires a lot of simulations of rather long regions to obtain enough simulated TEs, and the choice of relevant summary statistics is non-trivial (although abcrf can handle many summary statistics). There is also the issue of scale: one can simulate a fraction of the genome, retrieve summary statistics for the corresponding subset of simulated TEs, and infer on observed TEs drawn from the entire genome ; but there is a risk of distortion between simulations and observations: a subset of the genome is not always enough to cover the entire variance in recombination rates or gene content, assuming that we want to incorporate this in our simulations. We can of course pick carefully regions that are representative, or infer on genomic subsets with a similar distribution of recombination rate and gene densities as the windows used for training, but this may become quickly cumbersome and unpractical. The development of new deep learning methods may provide a more elegant and efficient way to retrieve parameters while learning to generalize to new configurations and number of observations in the observed data, while keeping the number of simulations relatively low (in the thousands range).
Deep learning: motivation and requirements
Two recent papers recently drew my attention : popgen-npe (Min et al. 2026) and an attention-based approach of Sanchez and colleagues (Sanchez et al. 2026). Both aim to make deep learning more accessible and more readily generalisable to genome-wide data.
An ideal pipeline would:
- take a (phased or unphased, polarised or not) VCF of TEs and SNPs
- extract features relevant to estimate TE age, frequency, genotypes, and where available distance to genes or other functional features and local recombination rate
- compare these features with simulations spanning a broad range of selection coefficients, proportions of strongly versus mildly deleterious insertions, transposition bursts, and insertion preferences and selection that varies along the genome. An advantage of the
popgen-npeimplementation is the ability to obtain an estimate of uncertainties around parameter estimation using neural posterior estimates. I wanted to incorporate this in the pipeline. - be as computationally efficient as possible. Whole-genome SLiM simulations are slow even when only deleterious markers and a few thousand segregating sites are modelled, so the method should learn from a representative subset of the genome and generalise to the whole. Of particular interest in the Sanchez et al. attention-based paper was the ability for the algorithm to weigh demographic inference obtained from distinct different genomic windows to produce a consensus demographic history through a voting scheme. In the preprint, they train the model on 200 x 2Mb windows but infer using the entire human genome (nearly ten times more windows). It seemed like a valuable way of using on TEs drawn from a genome of 1Gb a model trained on, for example, 20 x 5Mb windows (100Mb).
- optionally support iterative refinement in the spirit of ABC-MCMC or the sequential schemes used in tools such as DILS (Fraïsse et al. 2021): resampling from a broad initial posterior and refining, without overfitting
- remain inspectable. Deep learning is often called a black box, and coming from classical ABC I wanted a method whose behaviour I could interrogate: what the summary statistics look like, how the simulations behave, and how well they fit the observed data. A colleague of mine, François Sabot, told me about how some deep-learning methods designed to distinguish between images of huskies and other dog breeds were using the snow background as a distinguishing feature. Having the option to implement even a coarse ABC inference seemed useful to me from that perspective.
Working with AI
Large language models have become common and debated as research tools over the past few years. I had previously used Claude (Opus 4.7 to 5.0 here) for auxiliary tasks, mainly producing figures, improving my written English, and improving on Python, so I began a conversation to build my understanding of deep learning. This developed into what I later learned is called “vibe coding”: I give instructions, then collect, check and test the code. Python is fairly easy to reverse-engineer, so catching errors is easier than writing code from scratch.
What really made coding these trials much easier was the assistance AI provided for coding functions wrapping the simulators and the inference pipelines together. In my experience, this has always been the most difficult part of bioinformatics pipeline, converting the output of one (simulation) program into the input for another (inference) portion. It also really helped documenting the pipeline, and to make it adapted to a HPC environment.
On the more dangerous side, the AI agent is not particularly creative, and the speed at which it produces and edit code creates a perverse incentive to always complexify the method before properly testing what is already there. It is also biased towards more common languages, which is fine since Python is the main language used for the deep learning part. However, the SLiM code needed several corrections, especially since some functions in Eidos look like, but are not exactly the same, as their R equivalent, which generated a few issues during the first trials (for example, setting the gamma distribution in SLiM/Eidos involves specifying the shape and the mean, while in R one specifies the shape and the rate). Using summary statistics and ABC approaches was also useful to identify bugs (produced by the AI) in the scripts: for example, the benchmarking simulator only sampled windows with lower recombination rates in an earlier version. Interestingly, the benchmark was still decent, although a PCA on summary statistics clearly revealed the drift of benchmark simulations compared to the ones used for training.
What is shown below is the result of roughly three months of work (first prompts beginning of May 2026) to obtain a functional pipeline, starting by building on the popgen-npe architecture (Min et al. 2026) (itself a wrapper around the excellent sbi (Tejero-Cantero et al. 2020)) and on ideas from dnadna (Sanchez et al. 2026) about using attention to generalize inference to an arbitrary number of TEs beyond the range used in training. At this stage, I haven’t tested this framework on biological data, and all the parts of the code dealing with polarization, masking or phasing are yet to be tested. If you think this work may bring something to your research, or if you have any comment, feel free to reach out.
The hardest part was not, perhaps, the deep learning itself, for which there are a lot of excellent tutorials now. Reading the wiki pages for the sbi python library helped me a lot to understand the basic concepts around flows or embeddings. I could also count on a few lectures I could find. For example, I can only recommend Fumagalli, Perez and Jay’s lecture on deep learning here: https://github.com/mfumagalli/ImaGene/blob/master/Tutorials/Workshop/lecture.pdf
A typical simulation-based deep learning pipeline is organized this way (again, popgen-npe is really useful to understand the process, which will be familiar to ABC practitioners):
A simulator produces training data. I use SLiM v4.3 (Haller and Messer 2023)here, and intend to upgrade to v5. I also use msprime (Baumdicker et al. 2022) as a fast, rough simulator whose results can be recalibrated using a few SLiM simulations (see below).
A processor extracts the relevant features from the simulations, be it an image of the alignment (CNN), information about local recombination rate, gene density, or frequency of TEs.
A network architecture/flow can be then specified, that will take these features as input, and can be trained to transform the initial features into a set that can be used to run either regression (what we do here) or classification tasks
Once trained, we can benchmark on an independent set of simulations. In our case, we benchmark over a range of TEs or 5 Mb genomic windows.
If the algorithm seems to perform reasonably well, we can then infer based on actual data.
Overall, the modelling decisions, the diagnostics, and the interpretation are mine, AI accelerated implementation and helped me reason through alternatives in a domain that I did not master. I typically worked from my preliminary knowledge of ABC to identify equivalents in deep learning. For example, by asking whether there was an equivalent for improving posterior estimates through an enrichment of simulations around observed summary statistics (like in ABC-MCMC), the agent identified Sequential Neural Posterior estimate as a rough equivalent.
From images to sets of TEs.
The initial design used images of the TE genotype matrix together with flanking regions, the CNN model below. This is close to the model in popgen-npe, which in turn follows the exchangeable CNN model in Chan et al (Chan et al. 2018). I first chose this mostly because it seemed elegant to have an algorithm that could directly extract relevant information and associations between genomic features and TEs from the alignment.
Its drawbacks are essentially practical. Training is slow and memory-hungry, and the posteriors it returns are rather broad (a behaviour I still need to investigate). The tests below also use small effective population sizes. Production settings would be larger, but even then, modelling Drosophila populations (Ne on the order of two million) directly in SLiM is prohibitive and would require rescaling, reducing effective sizes and times while inflating mutation, insertion, and selection rates by the same factor. There is an excellent literature on doing this correctly (Marsh, Kaushik, and Johri 2026), which cautions against using windows larger than ~50-200 kb under two-orders-of-magnitude rescaling (it actually depends on the number of crossing-over events per simulated sequence/genome, see their Figure 1). In that regime the CNN approach, which needs a reasonably long genomic window, may be poorly suited.
The alternative, which is also more amenable to ABC, is to compute summary statistics on individual TEs, obtainable from much shorter, independent windows, and treat them as an unordered set of observations. We can however keep simulating genomic blocks, which may help keeping Linkage Disequilibrium between TEs, in the same way one can do block-bootstrap in population genetic inference to preserve correlations in the data. Here, I kept the same design throughout, with 5Mb windows, and sampled TEs from these windows, but the approach would also work when simulating many shorter windows.
The current pipeline extracts, for the SNPs flanking each TE, proxies for haplotype extension, local Tajima’s D for the carrier and non-carrier haplotype groups, local heterozygosity, and related quantities. These can be aggregated for ABC (the whole population of TEs is summarised: mean and quantiles of the recombination rate at insertion sites, of frequency, of nearby heterozygosity, and so on), which also mitigates the black-box problem, since I can inspect what the simulations produce and how well they match the data. The same per-TE summaries can equally be fed to a multi-layer perceptron (MLP).
A key point is that we never simulate an entire genome. The attention-based work by Sanchez et al. suggests a way to improve inference as more TEs are added without simulating their full genomic context. The design I settled on applies attention gates to the features produced by the CNN or the MLP, trains on a variable number of TEs so that differences in variance are less of an issue later on during inference, and extrapolates at inference to the observed data. This appears to work (see below): adding more TEs, even beyond the number of windows used in training continues to improve inference, most noticeably for the selection coefficient.
Technically this is not a transformer: the gates are computed from each TE’s own features, with no attention between TEs. It belongs to the Deep Sets family (Zaheer et al. 2018), using gated attention pooling in the sense of attention-based multiple-instance learning. Each TE contributes independently, adding more does not perturb those already there.
There are three main strategies that I tried, which all give roughly similar results, though there seems to be a slight advantage in using softmax attention gates, compared to sigmoid attention, or deepsets, which does not weigh individual TEs, but examines the distribution of per-TE features before feeding to the trunk. I also tried (mostly for myself) to provide an explanation to all this, which can be found as an annex.
Architectures
The pipeline implements six models, crossing three input/encoder strategies with two parameter sets. The two parameter sets are the DFE model, a gamma distribution of fitness effects (shape and mean \(|s|\)), the burst parameters (timing, duration, amplitude), and the global insertion rate ; and the context model, which additionally estimates how selection and insertion vary with local recombination through two slopes, \(\alpha_r\) (selection) and \(\gamma_r\) (insertion).
The parameterisation of the recombination dependence is
\[ s_{\text{local}}(p) = \exp\!\big(\log s_{\text{global}} + \alpha_r \, \phi_r(p)\big), \qquad \lambda_{\text{local}}(p) = \lambda_{\text{global}} \, \exp\!\big(\gamma_r \, \phi_r(p)\big), \]
where \(\phi_r(p)\) is a normalised log local recombination rate. Note that \(s_{\text{global}}\) and \(s_{\text{local}}\) are positive in the equation, but are then turned negative (purifying selection).
In the context model, the \(\lambda_{\text{global}}\) that we report in the tests below is the insertion rate expected for a specific median value of local recombination rate (RECOMB_TYPICAL in the config_base.py file), not the average insertion rate for the whole genome.
Note that we could very well use linear or quadratic functions, or use another feature (maybe chromatin accessibility, or gene density). I use recombination since it is usually correlated to gene density and chromatin accessibility, and is in itself a likely determinant of TE fitness due to ectopic recombination (Kent, Uzunović, and Wright 2017). Editing this should not require too much work, and most of the pipeline is already in place since I include distance to nearest gene as a feature used by the algorithm.
We fix the burst duration over 100 generations before rescaling. This is mostly because a 100 generations burst with a 100X increase from the baseline rate is not easy to distinguish from a 500 generations burst with a 20X increase at this scale, and I wanted to start testing relatively simple models. Note however that this is also a parameter that can be estimated in the current pipeline.
The three strategies are: (1) the CNN on genotype-matrix images with flanks; (2) normalised summary statistics aggregated across all TEs, fed to ABC; and (3) per-TE summary statistics treated as an unordered set, encoded by a shared MLP with attention gates. All results below use the same attention scheme (softmax gates, see Annex 2), which seems to produce slightly better results, but differences between the three (softmax, sigmoid, deepsets) are rather small.
(a) CNN + attention + DFE
(b) Normalised sumstats + ABC / neuralnet + DFE
flowchart TD
A[Per-TE summary statistics] --> B[Aggregate ACROSS TEs<br/>mean, sd, quantiles per feature]
B --> C{ABC-RF /-neuralnet / rejection}
C --> F[DFE parameters<br/>gamma shape, mean s,<br/>burst, insertion rate]
(c) Per-TE sumstats + MLP + attention + DFE
flowchart TD A[Per-TE summary statistics<br/>unordered set] --> B[Shared per-TE MLP encoder] B --> C[attention gates/deepsets + quantile pooling<br/>variable number of TEs K] C --> D[Neural posterior estimator<br/>normalising flow, sbi] D --> E[DFE parameters<br/>gamma shape, mean s,<br/>burst, insertion rate]
(d) CNN + attention + context
flowchart TD A[TE genotype matrix + flanking SNPs<br/>+ local recombination map] --> B[Two-branch CNN<br/>state + continuous, recomb-aware] B --> C[Per-window embeddings] C --> D[attention pooling/deepsets<br/>variable number of TEs K] D --> E[Neural posterior estimator<br/>normalising flow, sbi] E --> F[Context parameters<br/>DFE + burst + insertion rate,<br/>plus alpha_r selection slope<br/>and gamma_r insertion slope]
(e) Normalised sumstats + ABC / neuralnet + context
flowchart TD
A[Per-TE summary statistics<br/>incl. local recombination] --> B[Aggregate ACROSS TEs<br/>mean, sd, quantiles; recomb-stratified]
B --> C{ABC-RF /-neuralnet / rejection}
C--> D[Context parameters<br/>+ alpha_r, gamma_r]
(f) Per-TE sumstats + MLP + attention + context
flowchart TD A[Per-TE summary statistics<br/>unordered set, incl. local recombination] --> B[Shared per-TE MLP encoder] B --> C[Attention gates/deepsets + quantile pooling<br/>+ recombination summary features] C --> D[Neural posterior estimator<br/>normalising flow, sbi] D --> E[Context parameters<br/>DFE + burst + insertion rate,<br/>plus alpha_r and gamma_r]
Simulation set-up
The preliminary tests use a population of constant effective size with background selection at coding sites, driven by a recombination and gene map from a pilot project on the argan tree (Sideroxylon spinosum). The map was obtained using pyrho (Spence and Song 2019). The intent is a realistic picture of how gene density and recombination co-vary. This remains exploratory, and I will move to more complex/realistic demographic models later. To speed up the 5000 forward simulations that are used here, all quantities are rescaled by a factor of 10: effective sizes and times are divided by ten, and mutation, recombination, insertion, and selection rates are multiplied by ten. The parameter combinations are summarised below.
| Parameter | Symbol | Unscaled | Rescaled (10×) |
|---|---|---|---|
| Effective population size | \(N_e\) | 5000 diploids | 500 diploids |
| SNP mutation rate | \(\mu\) | \(1.25\times10^{-8}\) | \(1.25\times10^{-7}\) |
| Simulated duration | \(T\) | 50 000 gen | 5000 gen (\(=10\,N_e\)) |
| Window size | \(L\) | 5 Mb | — |
| Windows per replicate | — | 20 | — |
| Recombination / gene map | — | argan (S. spinosum). Recombination on average of the same order of magnitude as mutation rate. | — |
| Background sel. shape | BG_DEL_SHAPE |
0.3 | 0.3 |
| Background sel. mean \(s\) | BG_DEL_MEAN_S |
0.01 | 0.1 |
| Background sel. mutation rate | BG_DEL_MU |
\(0.85\times10^{-7}\) | \(0.85\times10^{-6}\) |
At the moment, tests are conducted on a set of 50 diploid individuals (100 haplotypes). Future tests will explore the robustness of methods to low sample size and lack of phase.
Under these background-selection settings roughly two-thirds of new coding mutations are deleterious. Simulations run for ten rescaled \(N_e\) generations to approach equilibrium. This is of course extremely simplified, a proper benchmark would use more realistic values for \(N_e\) and test various demographies.
The two recombination-dependence slopes of the context model, \(\alpha_r\) and \(\gamma_r\), may use a visual intuition. The figure below shows how the local selection coefficient and the local insertion intensity vary across four orders of magnitude of recombination, for two extreme (opposite-sign) values of each slope.
Code
#| label: fig-recomb-scaling
#| fig-cap: "Local selection coefficient (left) and local insertion intensity (right) as a function of the local recombination rate, for two extreme slope values. alpha_r governs the selection panel and gamma_r the insertion panel. A positive slope concentrates the effect in high-recombination regions; a negative slope does the opposite. phi_r is the normalised log local recombination rate, clipped to [-2, 2]."
#| warning: false
#| message: false
#| fig-width: 9
#| fig-height: 4
library(ggplot2)
library(dplyr)
library(tidyr)
library(knitr)
knitr::opts_chunk$set(
cache = FALSE
)
## Normalised recombination feature, matching the pipeline:
## phi_r(r) = clip( log10(r / RECOMB_TYPICAL) / RECOMB_LOG_SCALE, -2, 2 )
RECOMB_TYPICAL <- 1e-7 # "typical" recombination rate -> phi_r = 0
RECOMB_LOG_SCALE <- 1 # one log10 unit -> one unit of phi_r
feature_r <- function(r) pmin(pmax(log10(r / RECOMB_TYPICAL) / RECOMB_LOG_SCALE, -2), 2)
## Global baselines (illustrative)
log_s_global <- log(0.01) # mean selection coefficient s = 0.01
lambda_global <- 1 # baseline insertion intensity (relative units)
## Recombination axis and two extreme slopes
r <- 10^seq(-9, -5, length.out = 300)
extremes <- c(-2, 2)
df <- expand.grid(r = r, slope = extremes) |>
mutate(
phi = feature_r(r),
Selection = exp(log_s_global + slope * phi), # s_local (alpha_r = slope)
Insertion = lambda_global * exp(slope * phi) # lambda_local (gamma_r = slope)
) |>
pivot_longer(c(Selection, Insertion),
names_to = "quantity", values_to = "value") |>
mutate(
quantity = factor(quantity, levels = c("Selection", "Insertion")),
slope = factor(slope, labels = c("slope = -1", "slope = +1"))
)
ggplot(df, aes(r, value, colour = slope, linetype = slope)) +
geom_line(linewidth = 1) +
facet_wrap(~ quantity, scales = "free_y",
labeller = as_labeller(c(Selection = "Selection (alpha_r)",
Insertion = "Insertion (gamma_r)"))) +
scale_x_log10() +
scale_y_log10() +
scale_colour_manual(values = c("#B3452F", "#378ADD")) +
labs(x = "Local recombination rate (per bp per generation)",
y = "Local value (log scale)",
colour = NULL, linetype = NULL) +
theme_minimal(base_size = 12) +
theme(legend.position = "bottom")
We can also show where along our argan tree genome TEs would fall or be counterselected. For the context model, this is something useful to check.
Code
# Model (matching te_dfe/config_base.py):
# phi_r(p) = clip( log10(r(p) / RECOMB_TYPICAL) / RECOMB_LOG_SCALE, -2, 2 )
# lambda_local = lambda_global * exp( gamma_r * phi_r ) # insertion intensity
# s_local = exp( log_s_global + alpha_r * phi_r ) # selection against TEs
#
#Colour is on a log scale because both quantities span ~3000-fold across the
# map at |slope| = 2 (exp(+/-4)).
library(ggplot2)
library(dplyr)
library(tidyr)
## ---- CONFIG: mirror te_dfe/config_base.py -----------------------------------
RECOMB_MAP_PATH <- "recomb_map.txt" # HapMap: Chromosome Position(bp) Rate(cM/Mb) Map(cM)
CHROMS <- NULL # NULL = all chromosomes in the file
RECOMB_TYPICAL <- 1e-7 # r giving phi_r = 0
RECOMB_LOG_SCALE <- 1 # one log10 unit of r = one unit of phi_r
LOG_S_GLOBAL <- log(0.01) # baseline selection coefficient (s = 0.01)
LAMBDA_GLOBAL <- 1 # baseline insertion intensity (relative units)
BIN_BP <- 1e5 # colour-bar resolution
GAMMA_VALUES <- c(-2, 2) # extreme insertion slopes (Figure 1)
ALPHA_VALUES <- c(-2, 2) # extreme selection slopes (Figure 2)
feature_r <- function(r) {
pmin(pmax(log10(r / RECOMB_TYPICAL) / RECOMB_LOG_SCALE, -2), 2)
}
## ---- Read the HapMap map ----------------------------------------------------
## Rate is cM/Mb; Morgan/bp = rate * 1e-8 (same conversion as
## te_dfe/demography.py::load_hapmap_recomb).
raw <- read.table(RECOMB_MAP_PATH, header = FALSE, fill = TRUE,
stringsAsFactors = FALSE,
col.names = c("chrom", "pos", "rate_cM_Mb", "map_cM"))
raw <- raw[!is.na(suppressWarnings(as.numeric(raw$pos))), ] # drop header row
raw$pos <- as.numeric(raw$pos)
raw$rate_cM_Mb <- as.numeric(raw$rate_cM_Mb)
raw$chrom <- as.character(raw$chrom)
if (!is.null(CHROMS)) raw <- raw[raw$chrom %in% CHROMS, ]
## Natural order: chr2 before chr10 (a plain sort would invert them).
uc <- unique(raw$chrom)
chrom_num <- suppressWarnings(as.numeric(gsub("[^0-9]", "", uc)))
chrom_lev <- uc[order(ifelse(is.na(chrom_num), Inf, chrom_num), uc)]
## Step segments, computed within each chromosome so none straddles a boundary.
bins <- raw |>
mutate(chrom = factor(chrom, levels = chrom_lev)) |>
arrange(chrom, pos) |>
group_by(chrom) |>
mutate(seg_start = lag(pos, default = 0),
seg_end = pos,
seg_len = seg_end - seg_start,
r = rate_cM_Mb * 1e-8) |> # cM/Mb -> Morgan/bp
ungroup() |>
filter(seg_len > 0, is.finite(r)) |>
## bin to BIN_BP: length-weighted mean rate per bin
mutate(bin = floor(((seg_start + seg_end) / 2) / BIN_BP)) |>
group_by(chrom, bin) |>
summarise(r = sum(r * seg_len) / sum(seg_len), .groups = "drop") |>
mutate(mid_mb = (bin + 0.5) * BIN_BP / 1e6,
phi = feature_r(r))
cat(sprintf("map: %d chromosomes, %d bins, phi_r in [%.2f, %.2f]\n",
nlevels(droplevels(bins$chrom)), nrow(bins),
min(bins$phi), max(bins$phi)))map: 11 chromosomes, 5496 bins, phi_r in [-2.00, 0.92]Code
## ---- Shared plot skeleton ---------------------------------------------------
chrom_bars <- function(df, fill_name, title, subtitle, option = "viridis") {
ggplot(df, aes(x = mid_mb, y = chrom, fill = value)) +
geom_tile(height = 0.7, width = BIN_BP / 1e6) +
facet_wrap(~ panel, ncol = 1) +
scale_fill_viridis_c(trans = "log10", option = option, name = fill_name) +
scale_y_discrete(limits = rev(levels(df$chrom))) +
scale_x_continuous(expand = c(0.01, 0)) +
labs(x = "Position (Mb)", y = NULL, title = title, subtitle = subtitle) +
theme_minimal(base_size = 11) +
theme(panel.grid = element_blank(),
legend.position = "bottom",
legend.key.width = unit(2.5, "cm"),
strip.text = element_text(face = "bold"))
}
## ---- Figure 1: insertion intensity (gamma_r) --------------------------------
fig1_df <- tidyr::crossing(bins, gamma_r = GAMMA_VALUES) |>
mutate(value = LAMBDA_GLOBAL * exp(gamma_r * phi),
panel = paste0("gamma_r = ", sprintf("%+g", gamma_r)))
fig1 <- chrom_bars(
fig1_df,
fill_name = "Insertion intensity (relative to baseline)",
title = "Figure 1 — Where TEs insert preferentially",
subtitle = "lambda_local = lambda_global * exp(gamma_r * phi_r). Bright = TEs inserted more often.",
option = "viridis"
)
print(fig1)
Code
## ---- Figure 2: selection coefficient (alpha_r) ------------------------------
fig2_df <- tidyr::crossing(bins, alpha_r = ALPHA_VALUES) |>
mutate(value = exp(LOG_S_GLOBAL + alpha_r * phi),
panel = paste0("alpha_r = ", sprintf("%+g", alpha_r)))
fig2 <- chrom_bars(
fig2_df,
fill_name = "Selection coefficient s against TEs",
title = "Figure 2 — Where TEs are counter-selected",
subtitle = "s_local = exp(log_s_global + alpha_r * phi_r). Bright = selection removes TEs more efficiently.",
option = "magma"
)
print(fig2)
Code
## ggsave("fig1_insertion_gamma_r.png", fig1, width = 10, height = 6, dpi = 150)
## ggsave("fig2_selection_alpha_r.png", fig2, width = 10, height = 6, dpi = 150)Arguably, a model with strong selection against TEs in regions of lower recombination is less likely, but I was interested in checking whether \(\alpha_r\) and \(\gamma_r\) were indeed distinguishable.
Comparing the approaches
For each of the two models I report recovery of the target parameters against their true values, and I split between two sets of simulations: those with an intense transposition burst (> 100X the baseline), and a lower-intensity burst (< 100X). This is mostly because bursts that are not intense and old are difficult to distinguish from a flat influx of TEs. Note that we are using a rather low number of simulations, particularly for ABC inference (5000 simulations of 20 x 5Mb windows for training).
We can extract summary statistics about Tajima’s D, heterozygosity quantiles, haplotype extension for each simulation. The full list of sumstats is in annex. There are many summary statistics, but we can simplify them into a PCA and run that into a classical ABC framework. This can be a useful comparison to more advanced deep learning techniques, and makes it easier to explore the behaviour of our statistics depending on the parameters. It is also very helpful to spot bugs in the model, or determine whether we approximate our observed data well.
Training was done by randomly subsampling between 20 and 1000 TEs (capped at 1000 TEs collected over 20 windows), to expose the algorithm to varying amounts of variance. Of course, some statistics still reflect the original number of TEs in the dataset (for example, higher transposition rates give more TEs and therefore higher TE density values in flanking windows). We benchmark by running the trained algorithm on a set of 200 independent simulations of 200 x 5Mb windows that were not used for training. This allows testing the robustness of our inference pipelines to a number of TEs that is never encountered during training (ten times more). Note that 10,000 polymorphic TEs may be on the high side for any given lineage (let’s say ATHILA in plants, or any consensus of interest), although such numbers may be reached by pooling TEs into broader families (for example all Ty1/Copia together, although this may not make much sense biologically in some cases). The point is that the method is trained on 1000, which may be a more realistic number for (abundant) lineages in most genomes.
(a) Normalised sumstats + ABC/network + DFE
We can obtain the distribution of per-TE summary statistics using this python script, which runs on the datasets generated for each of the SLiM simulations (see also the git repository and the annex for more commands and explanations on the simulation + training commands):
Code
python 07_export_sumstat_summaries.py --model simple --dataset <path_to_collated_training_dataset> --out summaries_simple.csvWe also obtain the same summary statistics automatically for each 200 benchmarking datasets (themselves containing 200 x 5 Mb windows spanning a range of gene density and recombination rates) when we run the benchmarking scripts after training the deep learning algorithm. So we can now use ABC to examine our summary stats and compare.
Code
# ABC benchmark on the per-TE summary statistics.
library(abc)
library(HDInterval)
library(dplyr)
library(ggplot2)
library(patchwork)
library(readr); library(dplyr); library(stringr)
library(FactoMineR); library(data.table)
base="Results_final_blogpost/DFE_model_sumstats/"
ABC_METHOD <- "neuralnet" # "rejection" is ~10x faster for a first pass but not very efficient
ABC_TOL <- 0.05
CRED <- 0.90 # can be set to any interval deemed interesting, can be a way to check different slices of the posterior distribution. More of a sanity check.
NPCS <- 10
All <- na.omit(read.csv(file.path(base,"../summaries_simple.csv"), header = TRUE, fill = TRUE))
parameters <- All[, 1:6]
##mean= alpha/beta with beta the rate and alpha the shape. so rate is shape/mean
parameters$p_strong <- pgamma(100, shape = exp(parameters$log_gamma_shape), rate = exp(parameters$log_gamma_shape)/exp(parameters$log_mean_s)) - pgamma(1e-2, shape = exp(parameters$log_gamma_shape), rate = exp(parameters$log_gamma_shape)/exp(parameters$log_mean_s))
parameters$p_weak <- pgamma(1e-2, shape = exp(parameters$log_gamma_shape), rate = exp(parameters$log_gamma_shape)/exp(parameters$log_mean_s)) - pgamma(1e-3, shape = exp(parameters$log_gamma_shape), rate = exp(parameters$log_gamma_shape)/exp(parameters$log_mean_s))
parameters$p_nn <-pgamma(1e-3, shape = exp(parameters$log_gamma_shape), rate = exp(parameters$log_gamma_shape)/exp(parameters$log_mean_s)) - pgamma(0, shape = exp(parameters$log_gamma_shape), rate = exp(parameters$log_gamma_shape)/exp(parameters$log_mean_s))
train_df <- All[, 8:ncol(All)]
# drop zero-variance / constant columns and column about the number of TEs per simulation (which is clearly different here by one order of magnitude)
train_df <- train_df %>% dplyr::select(where(~ if (is.numeric(.)) var(., use = "complete.obs") != 0 else length(unique(.)) > 1))
res.pca <- PCA(train_df, graph = FALSE, ncp = 10)
output_coord <- data.frame(res.pca$ind$coord)
output_coord$Category <- ifelse(exp(parameters$log_mean_s) > 0.01,
"strong sel", "weaker sel")
pnames <- colnames(parameters)
P <- length(pnames)
train_pcs <- output_coord[, paste0("Dim.", seq_len(NPCS))]
# --- PCA scatter of training vs benchmark, with the CORRECT projection ---
observed_full <- read.csv(file.path(base, "benchmark_NTE10000_abc_summaries.csv"), header = TRUE)
observed_full <- na.omit(observed_full) # drop rows first -> alignment holds
obs_sumstat <- observed_full[, colnames(train_df), drop = FALSE] # same order
observed_df <- as.data.frame(predict(res.pca, newdata = obs_sumstat)$coord)
observed_df$Category <- "Observed"
print(
ggplot(rbind(output_coord, observed_df),
aes(Dim.1, Dim.2, color = Category)) +
geom_point(size = 3, alpha = 0.8) + theme_minimal() +
scale_color_brewer(palette = "Set1"))
The good thing is that the “observed” summary statistics fall within the cloud of training simulations, so using normalized statistics works reasonably well and is not too sensitive to the change in the number of TEs, although it would be more rigorous to compare simulations producing the same number (so 20 windows against 20 windows, not 200).
We can dig a bit more and examine the correlation of summary statistics with distinct parameters.
Code
library(dplyr); library(ggplot2)
param_cols <- colnames(All[, 1:6])
target <- "log_lambda" # parameter to screen
method <- "spearman" # "spearman" (robust) or "pearson"
partial <- TRUE # control for the other parameters?
top_n <- 50
stopifnot(target %in% names(All))
# stats = everything that isn't a parameter (n_valid_te kept as a candidate stat)
stat_cols <- setdiff(names(All), param_cols)
stat_cols <- stat_cols[sapply(All[stat_cols], is.numeric)]
v <- sapply(All[stat_cols], var, na.rm = TRUE) # drop constant / all-NA cols
stat_cols <- stat_cols[is.finite(v) & v > 0]
y <- All[[target]]
if (!partial) {
cors <- sapply(stat_cols, function(s)
suppressWarnings(cor(y, All[[s]], method = method, use = "pairwise.complete.obs")))
} else {
others <- setdiff(param_cols, target)
others <- others[sapply(All[others], function(x) var(x, na.rm = TRUE) > 0)] # drop fixed params
ctrl <- as.matrix(All[others])
resid_of <- function(z) residuals(lm(z ~ ctrl, na.action = na.exclude))
ry <- resid_of(y)
cors <- sapply(stat_cols, function(s)
suppressWarnings(cor(ry, resid_of(All[[s]]), method = method, use = "pairwise.complete.obs")))
}
corr_tbl <- tibble(stat = stat_cols, cor = as.numeric(cors)) |>
filter(is.finite(cor)) |> arrange(desc(abs(cor)))
print(corr_tbl, n = top_n)# A tibble: 129 × 2
stat cor
<chr> <dbl>
1 n_valid_te 0.707
2 local_te_density__mean 0.650
3 local_te_density__q50 0.648
4 local_te_density__q75 0.645
5 local_te_density__q25 0.644
6 local_te_density__q90 0.642
7 local_te_density__q10 0.639
8 local_te_density__max 0.636
9 local_te_density__sd -0.540
10 local_te_density__min 0.216
11 tajD_carrier__max 0.210
12 he_carrier_sd__max 0.204
13 tajD_carrier__min -0.192
14 he_noncarrier_sd__max 0.189
15 tajD_noncarrier__min -0.188
16 local_recomb__max 0.180
17 tajD_noncarrier__max 0.180
18 local_snp_density__min -0.179
19 he_noncarrier_mean__max 0.176
20 hapext_noncarrier__min -0.172
21 n_flank_snps_norm__min -0.168
22 local_snp_density__q90 -0.165
23 he_carrier_mean__max 0.163
24 local_snp_density__q75 -0.160
25 hapext_noncarrier__max 0.149
26 n_flank_snps_norm__max 0.142
27 local_snp_density__q50 -0.133
28 hapext_carrier__max 0.132
29 te_freq__max 0.126
30 he_noncarrier_sd__min -0.122
31 local_snp_density__mean -0.111
32 he_noncarrier_mean__min -0.109
33 n_flank_snps_norm__q75 -0.103
34 n_flank_snps_norm__mean -0.0932
35 n_flank_snps_norm__q50 -0.0839
36 n_flank_snps_norm__q25 -0.0790
37 local_snp_density__q25 -0.0787
38 n_flank_snps_norm__q10 -0.0749
39 hapext_carrier__q75 -0.0718
40 local_snp_density__q10 -0.0672
41 n_flank_snps_norm__q90 -0.0632
42 he_noncarrier_mean__sd 0.0607
43 hapext_carrier__q90 -0.0605
44 local_recomb__q25 0.0587
45 tajD_noncarrier__sd 0.0544
46 he_noncarrier_sd__sd 0.0535
47 hapext_carrier__mean -0.0524
48 he_noncarrier_mean__q10 -0.0504
49 tajD_noncarrier__q10 -0.0492
50 hapext_carrier__q50 -0.0483
# ℹ 79 more rowsCode
plot_tbl <- corr_tbl |> slice_head(n = top_n) |>
mutate(stat = factor(stat, levels = rev(stat)))
ggplot(plot_tbl, aes(cor, stat, color = cor > 0)) +
geom_segment(aes(x = 0, xend = cor, yend = stat), linewidth = 0.6) +
geom_point(size = 2.5) +
geom_vline(xintercept = 0, linewidth = 0.3) +
scale_color_manual(values = c(`TRUE`="#2c7fb8", `FALSE`="#de2d26"), guide = "none") +
labs(x = sprintf("%s %scorrelation with %s", method,
if (partial) "partial " else "", target),
y = NULL, title = sprintf("Stats most associated with %s", target)) +
theme_minimal(base_size = 11)