Quickstart

The easiest way to get started is by using the BaseInference class, which infers the DFE from a single pair of frequency spectra, one neutral and one selected. In this example we create Spectrum objects holding the SFS counts and pass them to BaseInference. The number of monomorphic sites must be specified: the first and last entries of the counts are the numbers of sites where the ancestral and derived allele is fixed, respectively. By default, only the deleterious part of the DFE is inferred (cf. fixed_params).

Python

import fastdfe as fd

inf = fd.BaseInference(
    sfs_neut=fd.Spectrum([177130, 997, 441, 228, 156, 117, 114, 83, 105, 109, 0]),
    sfs_sel=fd.Spectrum([797939, 1329, 499, 265, 162, 104, 117, 90, 94, 119, 0]),
    do_bootstrap=False
)

inf.run();

INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [--:--<--:--, --it/s]
INFO:Optimization: Optimizing 2 parameters: [all.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 3.9, all.b: 0.1305 ± 1.6e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 4.1e-09} (best_run ± std_across_runs)

R

library(fastdfe)
fd <- load_fastdfe()

inf <- fd$BaseInference(
  sfs_neut = fd$Spectrum(c(177130, 997, 441, 228, 156, 117, 114, 83, 105, 109, 0)),
  sfs_sel = fd$Spectrum(c(797939, 1329, 499, 265, 162, 104, 117, 90, 94, 119, 0)),
  do_bootstrap = FALSE
)

sfs_modelled <- inf$run()

INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [--:--<--:--, --it/s]
INFO:Optimization: Optimizing 2 parameters: [all.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 3.9, all.b: 0.1305 ± 1.6e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 4.1e-09} (best_run ± std_across_runs)

fastdfe uses maximum likelihood estimation (MLE) to find the DFE. By default, 10 local optimization runs are carried out to make sure a reasonably good global optimum has been found. The DFE furthermore needs to be parametrized, where GammaExpParametrization is used by default. The standard deviation across optimization runs is also reported to give an idea of the reliability of the estimates. In this case, the standard deviations are low, indicating that the estimates are stable.

We can now plot the inferred DFE in discretized form (cf. plot_discretized()).

Python

inf.plot_discretized();

R

p <- inf$plot_discretized()

We can also plot a comparison of the selected modelled and observed SFS (cf. plot_sfs_comparison()).

Python

inf.plot_sfs_comparison();

R

p <- inf$plot_sfs_comparison()

Bootstrapping

To quantify uncertainty, we can perform parametric bootstrapping (cf. bootstrap()).

Python

inf.bootstrap(n_samples=100)

inf.plot_discretized();

BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.021e+04 ± 4.1e+04, all.b: 0.1354 ± 0.021, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.74 ± 4.9, i_best_run: 0.46 ± 0.5, likelihoods_std: 0.2563 ± 2.3} (mean ± std)

R

bootstraps <- inf$bootstrap(n_samples = 100L)

p <- inf$plot_discretized()

BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.021e+04 ± 4.1e+04, all.b: 0.1354 ± 0.021, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.74 ± 4.9, i_best_run: 0.46 ± 0.5, likelihoods_std: 0.2563 ± 2.3} (mean ± std)

By default, 2 optimization runs are performed per bootstrap sample, taking the best result (cf. n_bootstrap_retries). The standard deviation across runs is computed for each bootstrap sample, and the average of these standard deviations across all samples is reported to summarize the uncertainty of the bootstrap estimates. In this case, the uncertainty is low, indicating reliable bootstrap estimates.