Quickstart

Contents

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. Note that we are required to specify the number of monomorphic sites (the last and first entries of the specified counts which correspond to the number of mono-allelic sites where the ancestral and derived allele is fixed, respectively). By default, only the deleterious part of the DFE is inferred (cf. fixed_params).

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
)

# run inference
inf.run();

INFO:Discretization: Precomputing linear DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:02<00:00,  3.20it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:01<00:00,  6.40it/s]
INFO:BaseInference: Successfully finished optimization after 18 iterations and 84 function evaluations, obtaining a log-likelihood of -35.44
INFO:BaseInference: Inferred parameters: {all.S_d: -3.389e+04, all.b: 0.1305, all.p_b: 0, all.S_b: 1, all.eps: 0}.
INFO:BaseInference: Standard deviations across runs: {all.S_d: 1.061e+04, all.b: 0.04187, all.p_b: 0, all.S_b: 0, all.eps: 0, likelihood: 7.282}.

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 bound. The DFE furthermore needs to parametrized where GammaExpParametrization is used by default. We also report the standard deviation across optimization runs 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()).

inf.plot_discretized();
../../_images/54e6b7d7882f87101007e6ddb0abc04e4ca42e06b95b1528d3bdd4811aca1623.png

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

inf.plot_sfs_comparison();
../../_images/eb42fa7b378cbf2382e801201411b43269e7fedcd35e90b5e962906672c306eb.png

Bootstrapping#

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

inf.bootstrap(n_samples=100)

# redo the plotting
inf.plot_discretized();

BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:18<00:00,  5.55it/s]
INFO:BaseInference: Standard deviations across bootstraps: {all.S_d: 3.862e+04, all.b: 0.02186, all.p_b: 0, all.S_b: 0, all.eps: 0, likelihood: 5.922, i_best_run: 0.5016, likelihoods_std: 0.04005}.
INFO:BaseInference: Average standard deviation of likelihoods across runs within bootstrap samples: 0.005547.
../../_images/1e5901f40ec902d16baf3ffcb451557874ee2f38bda3bfefdca6fb21ff349340.png

By default, we perform 2 optimization runs 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, we can see that the uncertainty is quite low indicating reliable bootstrap estimates.

Contents