DFE inference

Estimating the deleterious DFE

A short overview of basic DFE inference and bootstrapping is available in the quickstart guide. DFE inference requires one neutral and one selected SFS. In this example we use the bundled Betula pendula (silver birch) data. By default, bootstrapping is performed automatically, and the inference estimates only the deleterious component of the DFE using GammaExpParametrization.

import fastdfe as fd

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])

inf = fd.BaseInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel
)

# run inference
inf.run()

inf.plot_discretized();

INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  8.41it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 21.32it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 4.2, all.b: 0.1305 ± 1.7e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 3.2e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:06<00:00, 15.49it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.482e+04 ± 3.9e+04, all.b: 0.1331 ± 0.022, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.9 ± 5.9, i_best_run: 0.5 ± 0.5, likelihoods_std: 0.007777 ± 0.046} (mean ± std)

It is always good practice to check the variability of estimates across optimization runs to ensure stability. Here, both the standard deviations across initial runs and across runs within each bootstrap sample are low, indicating stable estimates. Individual runs and bootstrap results can be inspected in the corresponding dataframes (cf. runs and bootstraps).

inf.runs.select_dtypes('number')

	all.S_d	all.b	all.p_b	all.S_b	all.eps	all.h	likelihood
0	-33890.424207	0.130546	0.0	1.0	0.0	0.5	-35.437968
1	-33898.475526	0.130543	0.0	1.0	0.0	0.5	-35.437968
2	-33894.520478	0.130544	0.0	1.0	0.0	0.5	-35.437968
3	-33896.257305	0.130544	0.0	1.0	0.0	0.5	-35.437968
4	-33893.352524	0.130545	0.0	1.0	0.0	0.5	-35.437968
5	-33898.090546	0.130543	0.0	1.0	0.0	0.5	-35.437968
6	-33894.032960	0.130544	0.0	1.0	0.0	0.5	-35.437968
7	-33887.045393	0.130547	0.0	1.0	0.0	0.5	-35.437968
8	-33899.733825	0.130542	0.0	1.0	0.0	0.5	-35.437968
9	-33889.470185	0.130546	0.0	1.0	0.0	0.5	-35.437968

inf.bootstraps.select_dtypes('number')

	S_d	b	p_b	S_b	eps	h	likelihood	i_best_run	likelihoods_std	alpha
0	-100000.000000	0.115271	0.0	1.0	0.0	0.5	-48.666719	1	2.112557e-09	0.0
1	-23100.047602	0.132767	0.0	1.0	0.0	0.5	-39.454392	0	6.220944e-10	0.0
2	-30896.328327	0.126525	0.0	1.0	0.0	0.5	-44.116520	1	5.387903e-10	0.0
3	-100000.000000	0.117041	0.0	1.0	0.0	0.5	-40.332661	1	3.066069e-02	0.0
4	-27268.784916	0.134769	0.0	1.0	0.0	0.5	-38.642477	0	5.865388e-10	0.0
...	...	...	...	...	...	...	...	...	...	...
95	-59604.675726	0.121741	0.0	1.0	0.0	0.5	-38.589821	1	5.672405e-10	0.0
96	-3523.479918	0.172044	0.0	1.0	0.0	0.5	-33.468514	1	2.810907e-10	0.0
97	-100000.000000	0.114860	0.0	1.0	0.0	0.5	-44.168668	1	3.818051e-01	0.0
98	-5645.681207	0.160873	0.0	1.0	0.0	0.5	-34.451027	1	3.464038e-10	0.0
99	-12974.655531	0.146723	0.0	1.0	0.0	0.5	-36.749411	1	1.489667e-09	0.0

100 rows × 10 columns

We can also plot the parameter distributions across bootstrap samples to visualize uncertainty. We see that the mean strength of deleterious selection S_d often reaches the lower bounds of -1e5. Perhaps a different DFE parameterization or a more complex DFE model would be more appropriate here. It should also be noted that the spectra used in this example are far from exemplary, as they contain few SNPs and have a small sample size.

inf.bootstraps[['S_d', 'b']].hist(figsize=(6, 2));

We can also inspect how parameters covary.

inf.bootstraps.assign(S_d=inf.bootstraps.S_d.abs()).plot.hexbin('S_d', 'b', xscale='log', gridsize=15);

We observe a slight dependence between the mean S_d and the shape parameter b of GammaExpParametrization. This is because a large fraction of moderately deleterious mutations and a smaller fraction of strongly deleterious mutations can leave a similar signal in the SFS. Spectra with larger sample sizes might facilitate disentangling the two.

Estimating beneficial effects

Parameters can be held fixed during maximum-likelihood optimization, and this was already done internally in the example above. By default, fastDFE infers only the deleterious DFE, fixes the ancestral-allele misidentification rate eps, and assumes semi-dominant mutations (h = 0.5) (see fixed_params). Here, we estimate the full DFE, allowing for beneficial mutations by letting the parameters S_b and p_b of GammaExpParametrization vary, while eps and h remain fixed. The fixed parameters are wrapped in a dictionary under the key all, meaning these settings apply to all SFS types, which matters when running joint inference (cf. JointInference).

# create inference object
inf = fd.BaseInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel,
    fixed_params=dict(all=dict(eps=0, h=0.5))
)

# run inference
inf.run()

inf.plot_discretized();

INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  7.46it/s]
INFO:Optimization: Optimizing 4 parameters: [all.S_b, all.p_b, all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:01<00:00,  6.90it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {all.S_b: (0.0001, 100, 100)} and lower bound for {all.p_b: (0, 0.0043678716, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -9823 ± 4.3e+03, all.b: 0.1549 ± 0.13, all.p_b: 0.004368 ± 0.091, all.S_b: 100 ± 52, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -34.73 ± 0.24} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:20<00:00,  4.91it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.192e+04 ± 3.4e+04, all.b: 0.75 ± 2.1, all.p_b: 0.06205 ± 0.087, all.S_b: 58.15 ± 49, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -40.92 ± 5.8, i_best_run: 0.39 ± 0.49, likelihoods_std: 0.3645 ± 0.67} (mean ± std)

The inferred full DFE shows substantial uncertainty, which is expected with a small sample and few SNPs. This is most pronounced for the [-1, 0] and [0, 1] bins, in which mutations are effectively neutral and provide little signal. Adjusting the discretization intervals can help reveal the structure more clearly (cf. plot_discretized()).

inf.plot_discretized(intervals=[-np.inf, -100, -10, -1, 1, np.inf]);

Ancestral-allele misidentification

We can also adjust for ancestral-allele misidentification by letting parameter eps vary. eps is the probability that an allele is misidentified as derived when it is actually ancestral, and vice versa (cf. misidentify()). This can correct biases to the SFS caused by mis-polarization, but eps is somewhat difficult to interpret because it is applied simultaneously to both the neutral and selected SFS. In addition, eps assumes the fraction of ancestral misidentification to be constant across site classes, but in practise errors may differ across classes. Nevertheless, below, we infer the full DFE while allowing eps to vary.

inf = fd.BaseInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel,
    fixed_params=dict(all=dict(h=0.5))
)

inf.run();

INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  7.68it/s]
INFO:Optimization: Optimizing 5 parameters: [all.S_b, all.b, all.eps, all.S_d, all.p_b].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:01<00:00,  5.05it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.p_b: (0, 0, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1.004e+04 ± 7e+02, all.b: 0.1508 ± 0.0042, all.p_b: 0 ± 0.0065, all.S_b: 0.000167 ± 0.065, all.eps: 0.006854 ± 0.00012, all.h: 0.5 ± 0, likelihood: -34.63 ± 0.0061} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:28<00:00,  3.54it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.061e+04 ± 3.3e+04, all.b: 0.7611 ± 2.1, all.p_b: 0.07106 ± 0.09, all.S_b: 0.04142 ± 0.27, all.eps: 0.006649 ± 0.0068, all.h: 0.5 ± 0, likelihood: -40.64 ± 5.6, i_best_run: 0.53 ± 0.5, likelihoods_std: 0.01924 ± 0.074} (mean ± std)

inf.bootstraps[['eps']].hist(figsize=(4, 2));

eps is estimated to be rather low, indicating that ancestral-allele misidentification is not a major issue in this dataset, or, at least, that including it that does not significantly improve the model fit. We can also check this in a more principled way by performing a likelihood-ratio test as done below.

Nested model comparison

We can automatically check for the significance of include ancestral-allele misidentification and beneficial fitness affects using likelihood ratio tests. This is done with plot_nested_models(). The LRTs are performed by comparing the likelihood of the inferred DFE to the likelihood of a nested model where some parameters are held fixed. Alternatively, we can use compare_nested() to directly compare two nested models.

# set logging level to warning to avoid cluttering
fd.logger.setLevel('WARNING')

inf.plot_nested_models()

fd.logger.setLevel('INFO')

BaseInference>Performing inference: 100%|██████████| 10/10 [00:01<00:00,  6.02it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {all.S_b: (0.0001, 100, 100)} and lower bound for {all.p_b: (0, 0.0043680056, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
BaseInference>Performing inference: 100%|██████████| 10/10 [00:01<00:00,  5.24it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.p_b: (0, 0, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 25.52it/s]
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 16.20it/s]

Including ancestral allele misidentification or beneficial mutations does not appear to improve the fit significantly.

Dominance effects

By default, fastDFE assumes semi-dominance (h = 0.5), which is more or less appropriate depending on the organism and type of mutations considered. We can change the dominance coefficient to a different value of h if we believe this is more appropriate. However, in practise, h often depends on the strength of selection, with more deleterious mutations being more recessive. To model this, we can specify a callback function that returns the dominance coefficient as a function of the selection coefficient S = 4Ne.

In the example below, we use an exponential decay: h is about 0.4 for neutral mutations and approaches 0 for strongly deleterious ones. The callback also receives h itself, allowing the dominance function to be parametrized and optimized; for simplicity, this parameter is still called h. Its bounds can be set via bounds.

import numpy as np

inf = fd.BaseInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel,
    fixed_params=dict(all=dict(eps=0, h=0, p_b=0, S_b=1)),
    h_callback=lambda h, S: 0.4 * np.exp(-0.1 * abs(S))
)

inf.run();

INFO:Discretization: Precomputing DFE-SFS transformation for fixed dominance coefficients.
Discretization>Precomputing: 100%|██████████| 1809/1809 [00:32<00:00, 56.02it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 41.92it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
WARNING:BaseInference: The L1 residual comparing the modelled and observed SFS is rather large: `norm(sfs_modelled - sfs_observed, 1) / sfs_observed` = 0.157. This may indicate that the model does not fit the data well.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 0, all.b: 0.1401 ± 2.4e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0 ± 0, likelihood: -74.33 ± 1.3e-12} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:03<00:00, 29.31it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -1e+05 ± 0, all.b: 0.1396 ± 0.0032, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0 ± 0, likelihood: -83.91 ± 15, i_best_run: 0.55 ± 0.5, likelihoods_std: 2.238 ± 22} (mean ± std)

Let’s compare the inferred DFE under this dominance relationship to that of the default semi-dominant model.

inf2 = fd.BaseInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel,
    fixed_params=dict(all=dict(eps=0, h=0.5, p_b=0, S_b=1))
)

inf2.run();

INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  8.24it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 27.87it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 4.2, all.b: 0.1305 ± 1.7e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 3.2e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:05<00:00, 17.10it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.482e+04 ± 3.9e+04, all.b: 0.1331 ± 0.022, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.9 ± 5.9, i_best_run: 0.5 ± 0.5, likelihoods_std: 0.007777 ± 0.046} (mean ± std)

fd.DFE.plot_many([inf.get_dfe(), inf2.get_dfe()], labels=['partly recessive', 'h=0.5']);

We see that assuming mutations are partly recessive leads to a more deleterious inferred DFE since stronger selection is necessary to remove a similar amount of recessive mutations.

We can also let h vary when inferring the DFE (cf. the simulation guide).

Folded inference

To infer the DFE from a folded SFS, simply pass folded spectra to BaseInference. Folded inference is performed whenever the spectra are folded, i.e., when all entries where the derived allele is the major allele are zero. Folded spectra contain little information on beneficial mutations so we only infer the deleterious part of the DFE here.

import matplotlib.pyplot as plt
import numpy as np

# create inference object
inf = fd.BaseInference(
    sfs_neut=sfs_neut.fold(),
    sfs_sel=sfs_sel.fold()
)

# run inference
inf.run()

# plot the inferred DFE and the SFS comparison
ax1, ax2 = plt.subplots(ncols=2, figsize=(7, 3.2))[1]

inf.plot_discretized(ax=ax1, show=False, intervals=np.array([-np.inf, -100, -10, -1, 0]))
inf.plot_sfs_comparison(ax=ax2, show=False)

plt.tight_layout()
plt.show()

INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  8.07it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 25.12it/s]
INFO:BaseInference: Inference results: {all.S_d: -1.558e+04 ± 1.2, all.b: 0.1464 ± 1.4e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -21.58 ± 1.6e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [00:06<00:00, 14.61it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -3.681e+04 ± 3.9e+04, all.b: 0.1487 ± 0.027, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -25.07 ± 4.2, i_best_run: 0.51 ± 0.5, likelihoods_std: 0.001907 ± 0.013} (mean ± std)

Serialization

Inference objects can be serialized to JSON files for later use (cf. to_file()).

# save the inference object to a file
# we can unserialized the inference by using BaseInference.from_file
inf.to_file("out/serialized.json")

# we can also save a short summary to fa ile
inf.get_summary().to_file("out/summary.json")

Joint inference

fastDFE supports joint inference of several SFS types, where any parameters can be shared between types. In this example, we create a JointInference object with two types, where we share eps, the rate of ancestral misidentification and S_d, the mean selection coefficient for deleterious mutations (cf. GammaExpParametrization). For more complex stratifications, see the Parser) module.

# neutral SFS for two types
sfs_neut = fd.Spectra(dict(
    pendula=[177130, 997, 441, 228, 156, 117, 114, 83, 105, 109, 0],
    pubescens=[172528, 3612, 1359, 790, 584, 427, 325, 234, 166, 76, 31]
))

# selected SFS for two types
sfs_sel = fd.Spectra(dict(
    pendula=[797939, 1329, 499, 265, 162, 104, 117, 90, 94, 119, 0],
    pubescens=[791106, 5326, 1741, 1005, 756, 546, 416, 294, 177, 104, 41]
))

# create inference object
inf = fd.JointInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel,
    shared_params=[fd.SharedParams(types=["pendula", "pubescens"], params=["S_d"])]
)

# run inference
inf.run();

INFO:JointInference: Using shared parameters [SharedParams(params=['S_d'], types=['pendula', 'pubescens'])].
INFO:JointInference: Including covariates: {}.
INFO:JointInference: Running marginal inference for type 'all'.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  8.13it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 28.89it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 0, all.b: 0.1066 ± 2.2e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -45.09 ± 3.5e-12} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inferences for types ['pendula', 'pubescens'].
INFO:JointInference: Running marginal inference for type 'pendula'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 27.00it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 4.4, all.b: 0.1305 ± 1.8e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 3.3e-09} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'pubescens'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 10/10 [00:00<00:00, 28.84it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 6.4e+03, all.b: 0.1035 ± 0.00072, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -43.81 ± 0.094} (best_run ± std_across_runs)
INFO:JointInference: Running joint inference for types ['pendula', 'pubescens'].
INFO:Optimization: Optimizing 3 parameters: [pendula:pubescens.S_d, pubescens.b, pendula.b].
JointInference>Performing joint inference: 100%|██████████| 10/10 [00:01<00:00,  8.37it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {pendula:pubescens.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:JointInference: Inference results: {pendula.b: 0.1168 ± 0.00083, pendula.p_b: 0 ± 0, pendula.S_b: 1 ± 0, pendula.eps: 0 ± 0, pendula.h: 0.5 ± 0, pubescens.b: 0.1035 ± 0.00098, pubescens.p_b: 0 ± 0, pubescens.S_b: 1 ± 0, pubescens.eps: 0 ± 0, pubescens.h: 0.5 ± 0, pendula:pubescens.S_d: -1e+05 ± 8e+03, likelihood: -79.46 ± 0.12} (best_run ± std_across_runs)
INFO:JointInference: Bootstrapping type 'all'.
BaseInference>Bootstrapping 'all' (2 runs each): 100%|██████████| 100/100 [00:05<00:00, 19.88it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -8.665e+04 ± 2.7e+04, all.b: 0.1094 ± 0.0073, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -54 ± 7.5, i_best_run: 0.46 ± 0.5, likelihoods_std: 0.0002605 ± 0.0016} (mean ± std)
INFO:JointInference: Bootstrapping type 'pendula'.
BaseInference>Bootstrapping 'pendula' (2 runs each): 100%|██████████| 100/100 [00:05<00:00, 16.69it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.482e+04 ± 3.9e+04, all.b: 0.1331 ± 0.022, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.9 ± 5.9, i_best_run: 0.5 ± 0.5, likelihoods_std: 0.007777 ± 0.046} (mean ± std)
INFO:JointInference: Bootstrapping type 'pubescens'.
BaseInference>Bootstrapping 'pubescens' (2 runs each): 100%|██████████| 100/100 [00:04<00:00, 21.74it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -9.267e+04 ± 2.1e+04, all.b: 0.1049 ± 0.0055, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -53.15 ± 7.6, i_best_run: 0.4 ± 0.49, likelihoods_std: 3.296e-05 ± 0.00033} (mean ± std)
JointInference>Bootstrapping joint inference (2 runs each): 100%|██████████| 100/100 [00:16<00:00,  6.19it/s]
INFO:JointInference: Bootstrap summary: {pendula.b: 0.12 ± 0.0069, pendula.p_b: 0 ± 0, pendula.S_b: 1 ± 0, pendula.eps: 0 ± 0, pendula.h: 0.5 ± 0, pendula.S_d: -8.714e+04 ± 2.5e+04, pubescens.b: 0.1058 ± 0.0059, pubescens.p_b: 0 ± 0, pubescens.S_b: 1 ± 0, pubescens.eps: 0 ± 0, pubescens.h: 0.5 ± 0, pubescens.S_d: -8.714e+04 ± 2.5e+04, likelihood: -96.88 ± 9.2, i_best_run: 0.5 ± 0.5, likelihoods_std: 6.776 ± 68} (mean ± std)

JointInference both runs the joint inference and marginal inference where each type is inferred separately. To see this better we can plot the inferred parameters for the different inference types.

inf.plot_inferred_parameters();

marginal.pendula and marginal.pubescens are the marginal inferences for the respective type. marginal.all is the marginal inference obtaining by adding up the spectra of all types. joint.pendula and join.pubescens are the joint inferences for the respective type. We can see that eps and S_d are indeed shared between the two. The parameter alpha in the plot denotes the proportion of beneficial non-synonymous substitutions. Each marginal inference is a BaseInference object itself and can be accessed via inf.marginal_inferences[type].

We can now also investigate to what extent the inferred DFEs differ:

inf.plot_discretized();

Model comparison

We can obtain information about the goodness of fit achieved by sharing the parameter by performing a likelihood ratio test (cf. perform_lrt_shared()). This compares the likelihood of the joint inference with the product of the marginal likelihoods.

inf.perform_lrt_shared()

INFO:JointInference: Simple model likelihood: -79.46473650022017, Complex model likelihood: -79.24744560917105, Total degrees of freedom: 1, Parameters at boundary: 0.

0.5097492591276529

The test is not significant, indicating that the simpler model of sharing the parameters explains the data sufficiently well. Indeed, we do not observe a lot of differences between the inferred parameters of joint and the marginal inferences.

Covariates

JointInference also supports the inclusion of covariates associates with the different SFS types. This provides more powerful model testing and reduces the number of parameters that need to be estimated for the joint inference. For a more interesting example we stratify the SFS of B. pendula by the sites’ reference base as is described in more detail in the parser module.

# instantiate parser
p = fd.Parser(
    n=10,
    vcf="https://github.com/Sendrowski/fastDFE/"
        "blob/dev/resources/genome/betula/"
        "all.polarized.deg.subset.200000.vcf.gz?raw=true",
    stratifications=[fd.DegeneracyStratification(), fd.AncestralBaseStratification()]
)

# parse SFS
spectra: fd.Spectra = p.parse()

# visualize spectra
spectra.plot();

INFO:Parser: Using stratification: [neutral, selected].[A, C, G, T].
INFO:Parser: Loading VCF file
INFO:FileHandler: Using cached file at /var/folders/w6/18ktl5312413jw46btlxrh59fzxvt4/T/011b01ee5cec.all.polarized.deg.subset.200000.vcf.gz
INFO:FileHandler: Using cached file at /var/folders/w6/18ktl5312413jw46btlxrh59fzxvt4/T/011b01ee5cec.all.polarized.deg.subset.200000.vcf.gz
Parser>Counting sites: 200000it [00:02, 82178.86it/s]
Parser>Processing sites: 100%|██████████| 200000/200000 [03:30<00:00, 948.92it/s] 
INFO:PolyAllelicFiltration: Filtered out 154 sites.
INFO:DegeneracyStratification: Number of sites with valid type: 64083
INFO:AncestralBaseStratification: Number of sites with valid type: 64083
INFO:Parser: Skipped 983 sites without ancestral allele information.
INFO:Parser: Included 64083 out of 200000 sites in total from the VCF file.

We now create the inference object from the spectra. In this contrived example we make up some covariates that covary with S_d, the mean strength of negative selection. Covariates introduce a linear relationship by default but this can be modified by specifying a custom callback function (see Covariate).

# create inference object
inf = fd.JointInference(
    sfs_neut=spectra[['neutral.*']].merge_groups(1),
    sfs_sel=spectra[['selected.*']].merge_groups(1),
    covariates=[fd.Covariate(param='S_d', values=dict(A=1, C=2, T=3, G=4))],
    n_runs=50  # increase number of initial runs for stability
)

# run inference
inf.run();

INFO:JointInference: Parameters ['S_d'] have covariates and thus need to be shared. Adding them to shared parameters.
INFO:JointInference: Using shared parameters [SharedParams(params=['S_d'], types=['A', 'C', 'G', 'T'])].
INFO:JointInference: Including covariates: {'c0': {'param': 'S_d', 'values': {'A': 1, 'C': 2, 'T': 3, 'G': 4}}}.
INFO:JointInference: Running marginal inference for type 'all'.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  8.05it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 32.61it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 0, all.b: 0.1069 ± 2.5e-08, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -26.61 ± 1.5e-10} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inferences for types ['A', 'C', 'G', 'T'].
INFO:JointInference: Running marginal inference for type 'A'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 33.18it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
WARNING:BaseInference: The L1 residual comparing the modelled and observed SFS is rather large: `norm(sfs_modelled - sfs_observed, 1) / sfs_observed` = 0.200. This may indicate that the model does not fit the data well.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 0, all.b: 0.0857 ± 2e-08, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -22.72 ± 2.3e-11} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'C'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 25.80it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 3.1e+03, all.b: 0.1219 ± 0.0004, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -19.47 ± 0.00011} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'G'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 31.50it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 1.6e+04, all.b: 0.1221 ± 0.024, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -18.88 ± 0.43} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'T'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 32.98it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 8.8e+03, all.b: 0.1079 ± 0.0012, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -18.44 ± 0.0042} (best_run ± std_across_runs)
INFO:JointInference: Running joint inference for types ['A', 'C', 'G', 'T'].
INFO:Optimization: Optimizing 6 parameters: [A:C:G:T.S_d, T.b, C.b, A:C:G:T.c0, G.b, A.b].
JointInference>Performing joint inference: 100%|██████████| 50/50 [00:13<00:00,  3.69it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {A:C:G:T.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
WARNING:JointInference: Numerical optimization did not terminate normally, so the result might be unreliable. Consider adjusting the optimization parameters (increasing `gtol` or `n_runs`) or decreasing the number of optimized parameters.
INFO:JointInference: Inference results: {A.b: 0.0857 ± 3.3, A.p_b: 0 ± 0, A.S_b: 1 ± 0, A.eps: 0 ± 0, A.h: 0.5 ± 0, C.b: 0.1219 ± 1.9, C.p_b: 0 ± 0, C.S_b: 1 ± 0, C.eps: 0 ± 0, C.h: 0.5 ± 0, G.b: 0.1221 ± 0.52, G.p_b: 0 ± 0, G.S_b: 1 ± 0, G.eps: 0 ± 0, G.h: 0.5 ± 0, T.b: 0.1079 ± 2, T.p_b: 0 ± 0, T.S_b: 1 ± 0, T.eps: 0 ± 0, T.h: 0.5 ± 0, A:C:G:T.S_d: -1e+05 ± 4.9e+04, A:C:G:T.c0: 0 ± 5.2e+03, likelihood: -79.5 ± 4.1e+02} (best_run ± std_across_runs)
WARNING:BaseInference: The L1 residual comparing the modelled and observed SFS is rather large: `norm(sfs_modelled - sfs_observed, 1) / sfs_observed` = 0.200. This may indicate that the model does not fit the data well.
INFO:JointInference: Bootstrapping type 'all'.
BaseInference>Bootstrapping 'all' (2 runs each): 100%|██████████| 100/100 [00:04<00:00, 20.35it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -7.856e+04 ± 3.8e+04, all.b: 0.1168 ± 0.026, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.92 ± 6.4, i_best_run: 0.44 ± 0.5, likelihoods_std: 0.000985 ± 0.0099} (mean ± std)
INFO:JointInference: Bootstrapping type 'A'.
BaseInference>Bootstrapping 'A' (2 runs each): 100%|██████████| 100/100 [00:04<00:00, 22.51it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -7.758e+04 ± 4.1e+04, all.b: 0.1052 ± 0.044, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -91.89 ± 1.8e+02, i_best_run: 0.48 ± 0.5, likelihoods_std: 0.02762 ± 0.26} (mean ± std)
INFO:JointInference: Bootstrapping type 'C'.
BaseInference>Bootstrapping 'C' (2 runs each): 100%|██████████| 100/100 [00:04<00:00, 20.66it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.573e+04 ± 4.8e+04, all.b: 0.1812 ± 0.084, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -85.79 ± 1.5e+02, i_best_run: 0.47 ± 0.5, likelihoods_std: 0.004738 ± 0.027} (mean ± std)
INFO:JointInference: Bootstrapping type 'G'.
BaseInference>Bootstrapping 'G' (2 runs each): 100%|██████████| 100/100 [00:04<00:00, 20.75it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -6.275e+04 ± 4.7e+04, all.b: 0.1613 ± 0.066, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -109.6 ± 2.2e+02, i_best_run: 0.48 ± 0.5, likelihoods_std: 0.001265 ± 0.01} (mean ± std)
INFO:JointInference: Bootstrapping type 'T'.
BaseInference>Bootstrapping 'T' (2 runs each): 100%|██████████| 100/100 [00:05<00:00, 19.16it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.627e+04 ± 4.7e+04, all.b: 0.1616 ± 0.088, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -50.91 ± 82, i_best_run: 0.49 ± 0.5, likelihoods_std: 0.01103 ± 0.11} (mean ± std)
JointInference>Bootstrapping joint inference (2 runs each): 100%|██████████| 100/100 [00:43<00:00,  2.29it/s]
INFO:JointInference: Bootstrap summary: {A.b: 0.09628 ± 0.021, A.p_b: 0 ± 0, A.S_b: 1 ± 0, A.eps: 0 ± 0, A.h: 0.5 ± 0, A.S_d: -7.297e+04 ± 4e+04, A.c0: -235.5 ± 1.4e+03, C.b: 0.1376 ± 0.029, C.p_b: 0 ± 0, C.S_b: 1 ± 0, C.eps: 0 ± 0, C.h: 0.5 ± 0, C.S_d: -7.321e+04 ± 4e+04, C.c0: -235.5 ± 1.4e+03, G.b: 0.1386 ± 0.035, G.p_b: 0 ± 0, G.S_b: 1 ± 0, G.eps: 0 ± 0, G.h: 0.5 ± 0, G.S_d: -7.347e+04 ± 4e+04, G.c0: -235.5 ± 1.4e+03, T.b: 0.1226 ± 0.028, T.p_b: 0 ± 0, T.S_b: 1 ± 0, T.eps: 0 ± 0, T.h: 0.5 ± 0, T.S_d: -7.334e+04 ± 4e+04, T.c0: -235.5 ± 1.4e+03, likelihood: -506.9 ± 5.4e+02, i_best_run: 0.28 ± 0.45, likelihoods_std: 157 ± 2.1e+02} (mean ± std)

Let’s visualize the inferred parameters

inf.plot_inferred_parameters();

We observe that S_d shows little variation across the jointly inferred types, because it does not change linearly with respect to the arbitrary covariates specified. Indeed, the median of the covariate across all bootstrap replicates is close to zero. Note that covariates are named c0, c1, etc., by default.

inf.bootstraps['A.c0'].median()

0.0

Model comparison

We can perform a likelihood ratio test to see whether including the covariates produces a significantly better fit than simply sharing the parameter in question among the types (cf. perform_lrt_covariates()).

inf.perform_lrt_covariates()

INFO:JointInference: Using shared parameters [SharedParams(params=['S_d'], types=['A', 'C', 'G', 'T'])].
INFO:JointInference: Including covariates: {}.
INFO:JointInference: Running joint inference without covariates.
INFO:JointInference: Running marginal inference for type 'all'.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [00:01<00:00,  8.54it/s]
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 32.00it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 0, all.b: 0.1069 ± 2.5e-08, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -26.61 ± 1.5e-10} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inferences for types ['A', 'C', 'G', 'T'].
INFO:JointInference: Running marginal inference for type 'A'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 34.27it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
WARNING:BaseInference: The L1 residual comparing the modelled and observed SFS is rather large: `norm(sfs_modelled - sfs_observed, 1) / sfs_observed` = 0.200. This may indicate that the model does not fit the data well.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 0, all.b: 0.0857 ± 2e-08, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -22.72 ± 2.3e-11} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'C'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 29.24it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 3.1e+03, all.b: 0.1219 ± 0.0004, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -19.47 ± 0.00011} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'G'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 34.58it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 1.6e+04, all.b: 0.1221 ± 0.024, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -18.88 ± 0.43} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'T'.
INFO:Optimization: Optimizing 2 parameters: [all.b, all.S_d].
BaseInference>Performing inference: 100%|██████████| 50/50 [00:01<00:00, 36.12it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1e+05 ± 8.8e+03, all.b: 0.1079 ± 0.0012, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -18.44 ± 0.0042} (best_run ± std_across_runs)
INFO:JointInference: Running joint inference for types ['A', 'C', 'G', 'T'].
INFO:Optimization: Optimizing 5 parameters: [A:C:G:T.S_d, T.b, C.b, G.b, A.b].
JointInference>Performing joint inference: 100%|██████████| 50/50 [00:13<00:00,  3.66it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {A:C:G:T.S_d: (-100000, -100000, -0.01)} [(lower, value, upper)], but this might be nothing to worry about.
WARNING:JointInference: Numerical optimization did not terminate normally, so the result might be unreliable. Consider adjusting the optimization parameters (increasing `gtol` or `n_runs`) or decreasing the number of optimized parameters.
INFO:JointInference: Inference results: {A.b: 0.0857 ± 0.0009, A.p_b: 0 ± 0, A.S_b: 1 ± 0, A.eps: 0 ± 0, A.h: 0.5 ± 0, C.b: 0.1219 ± 0.0011, C.p_b: 0 ± 0, C.S_b: 1 ± 0, C.eps: 0 ± 0, C.h: 0.5 ± 0, G.b: 0.1221 ± 0.0012, G.p_b: 0 ± 0, G.S_b: 1 ± 0, G.eps: 0 ± 0, G.h: 0.5 ± 0, T.b: 0.1079 ± 0.00096, T.p_b: 0 ± 0, T.S_b: 1 ± 0, T.eps: 0 ± 0, T.h: 0.5 ± 0, A:C:G:T.S_d: -1e+05 ± 8e+03, likelihood: -79.5 ± 0.023} (best_run ± std_across_runs)
WARNING:BaseInference: The L1 residual comparing the modelled and observed SFS is rather large: `norm(sfs_modelled - sfs_observed, 1) / sfs_observed` = 0.200. This may indicate that the model does not fit the data well.
INFO:JointInference: Simple model likelihood: -79.50114061683993, Complex model likelihood: -79.50114061683993, Total degrees of freedom: 1, Parameters at boundary: 0.

1.0

As expected, the specified covariates do not improve the fit significantly.