DFE inference

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.

Python

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
)

inf.run()

inf.plot_discretized();

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

../_images/13300e83d3d44cbf1d8073e65749dccd17066635a08d483356dcc1eff9f4fe72.png

R

library(fastdfe)
fd <- load_fastdfe()

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

inf <- fd$BaseInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel
)

sfs_modelled <- inf$run()

p <- inf$plot_discretized()

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

../_images/09bc08a1c6229e3f576879d109cb4def9b658a73f22f3be6d50075d44e7a7dcf.png

It is 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).

Python

inf.runs.select_dtypes('number')

	all.S_d	all.b	all.S_b	all.h	likelihood
0	-33901.428952	0.130541	1.0	0.5	-35.437968
1	-33892.197631	0.130545	1.0	0.5	-35.437968
2	-33892.125695	0.130545	1.0	0.5	-35.437968
3	-33895.005810	0.130544	1.0	0.5	-35.437968
4	-33892.726790	0.130545	1.0	0.5	-35.437968
5	-33900.321830	0.130542	1.0	0.5	-35.437968
6	-33889.389000	0.130546	1.0	0.5	-35.437968
7	-33897.772583	0.130543	1.0	0.5	-35.437968
8	-33892.245338	0.130545	1.0	0.5	-35.437968
9	-33895.271683	0.130544	1.0	0.5	-35.437968

inf.bootstraps.select_dtypes('number').head(10)

	S_d	b	S_b	h	likelihood	i_best_run	likelihoods_std	omega
0	-100000.000000	0.114344	1.0	0.5	-41.253020	0	5.684342e-13	0.223217
1	-8855.243409	0.149063	1.0	0.5	-48.350324	0	3.788898e-10	0.211289
2	-10536.012806	0.146062	1.0	0.5	-38.581950	0	2.527258e-10	0.211919
3	-100000.000000	0.120063	1.0	0.5	-34.680724	1	5.684342e-14	0.208291
4	-100000.000000	0.118915	1.0	0.5	-58.803044	0	4.897894e-01	0.211200
5	-100000.000000	0.114051	1.0	0.5	-36.241569	0	8.526513e-14	0.224013
6	-3005.297995	0.176505	1.0	0.5	-46.800673	0	4.670540e-10	0.197781
7	-2255.046969	0.188293	1.0	0.5	-39.336834	0	3.521166e-10	0.189612
8	-3097.896134	0.182649	1.0	0.5	-42.141935	1	1.014939e-10	0.187047
9	-33572.101276	0.133442	1.0	0.5	-44.470665	1	2.493981e-05	0.205144

R

inf$runs[sapply(inf$runs, is.numeric)]

A data.frame: 10 × 7
all.S_d	all.b	all.p_b	all.S_b	all.eps	all.h	likelihood
<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
-33901.43	0.1305413	0	1	0	0.5	-35.43797
-33892.20	0.1305452	0	1	0	0.5	-35.43797
-33892.13	0.1305452	0	1	0	0.5	-35.43797
-33895.01	0.1305440	0	1	0	0.5	-35.43797
-33892.73	0.1305449	0	1	0	0.5	-35.43797
-33900.32	0.1305418	0	1	0	0.5	-35.43797
-33889.39	0.1305463	0	1	0	0.5	-35.43797
-33897.77	0.1305428	0	1	0	0.5	-35.43797
-33892.25	0.1305451	0	1	0	0.5	-35.43797
-33895.27	0.1305439	0	1	0	0.5	-35.43797

head(inf$bootstraps[sapply(inf$bootstraps, is.numeric)], 10)

A data.frame: 10 × 12
	S_d	b	p_b	S_b	eps	h	likelihood	i_best_run	likelihoods_std	alpha	omega	omega_a
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	-100000.000	0.1143444	0	1	0	0.5	-41.25302	0	5.684342e-13	0	0.2232168	0
2	-8855.243	0.1490629	0	1	0	0.5	-48.35032	0	3.788898e-10	0	0.2112888	0
3	-10536.013	0.1460621	0	1	0	0.5	-38.58195	0	2.527258e-10	0	0.2119193	0
4	-100000.000	0.1200634	0	1	0	0.5	-34.68072	1	5.684342e-14	0	0.2082915	0
5	-100000.000	0.1189155	0	1	0	0.5	-58.80304	0	4.897894e-01	0	0.2112004	0
6	-100000.000	0.1140510	0	1	0	0.5	-36.24157	0	8.526513e-14	0	0.2240126	0
7	-3005.298	0.1765053	0	1	0	0.5	-46.80067	0	4.670540e-10	0	0.1977813	0
8	-2255.047	0.1882931	0	1	0	0.5	-39.33683	0	3.521166e-10	0	0.1896122	0
9	-3097.896	0.1826490	0	1	0	0.5	-42.14193	1	1.014939e-10	0	0.1870472	0
10	-33572.101	0.1334418	0	1	0	0.5	-44.47066	1	2.493981e-05	0	0.2051444	0

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

Python

inf.bootstraps[['S_d', 'b']].hist(figsize=(8.5, 2.5), grid=False, xrot=30);

../_images/b95f544dc98786ef6e2d87e77c3a223eef90b270d4b2a81854d0fe63c4bd217f.png

R

par(mfrow = c(1, 2), mar = c(2.5, 3, 2, 1))
hist(inf$bootstraps$S_d, main = "S_d", xlab = "", col = "#1f77b4", border = "white")
hist(inf$bootstraps$b, main = "b", xlab = "", col = "#1f77b4", border = "white")

../_images/7223165c2fa4aa354866870e6e790e37a23371e7d4cc75dd21e34709785af578.png

We can also inspect how parameters covary.

Python

inf.bootstraps.assign(S_d=inf.bootstraps.S_d.abs()).plot.scatter('S_d', 'b', logx=True);

../_images/b2c2ef59579c16e14de5b6d494cdda4f79224dfaa0adafd7dfbb862784fc288e.png

R

par(mar = c(4, 4, 1, 1))
plot(abs(inf$bootstraps$S_d), inf$bootstraps$b, log = "x", xlab = "S_d", ylab = "b", pch = 16, col = "#1f77b4")

../_images/1314cb42821d7fc94d19ce28389daa61afa5fe1e3f78e9e668d594cfc61eda95.png

We observe a strong 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 grouped under the key all, meaning these settings apply to all SFS types, which matters when running joint inference (cf. JointInference).

Python

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

inf.run()

inf.plot_discretized();

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 4 parameters: [all.S_d, all.b, all.p_b, all.S_b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.0043678413, 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 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.365e+04 ± 3.6e+04, all.b: 0.6986 ± 2, all.p_b: 0.05515 ± 0.079, all.S_b: 52.34 ± 49, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -40.62 ± 4.5, i_best_run: 0.41 ± 0.49, likelihoods_std: 0.3548 ± 0.64} (mean ± std)

../_images/424d887110e336c6f7925fec955cde6d644d6b5ab8e2e9b44cedf9a6103344fd.png

R

inf <- fd$BaseInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel,
  fixed_params = list(all = list(eps = 0, h = 0.5))
)

sfs_modelled <- inf$run()

p <- inf$plot_discretized()

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 4 parameters: [all.S_d, all.b, all.p_b, all.S_b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.0043678413, 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 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.365e+04 ± 3.6e+04, all.b: 0.6986 ± 2, all.p_b: 0.05515 ± 0.079, all.S_b: 52.34 ± 49, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -40.62 ± 4.5, i_best_run: 0.41 ± 0.49, likelihoods_std: 0.3548 ± 0.64} (mean ± std)

../_images/4c77c0a86fdbbad5f327a6d0cfc2541c20da0ace714d046f153af7b2ecb962a1.png

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()).

Python

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

../_images/a8041d1749d91b93556bbb7cbd321fc740b4d6ad4a15ae00dcb3211b51a5178b.png

R

p <- inf$plot_discretized(intervals = c(-Inf, -100, -10, -1, 1, Inf))

../_images/1c657e689cd288f22725069c91afe2cb06bc6def63f8b39af05a24b0aab1055e.png

Divergence counts#

Besides polymorphism, fastdfe can incorporate divergence counts, the number of fixed differences (substitutions) to an outgroup, into the likelihood, much like polydfe. The last entry of an SFS is the fixed-derived (divergence) class. To make use of divergence, the divergence target sizes n_sites_div_neut and n_sites_div_sel are additionally passed to the inference. They are the numbers of mutational target sites over which divergence was counted, which may differ from the polymorphism target size.

When both spectra carry a separate divergence target size, divergence is included in the likelihood automatically (see include_divergence). The selected divergence then helps constrain the beneficial part of the DFE, and \(\alpha\), the proportion of beneficial substitutions, can be estimated McDonald–Kreitman style from the observed divergence.

Python

# the SFS runs from the monomorphic (ancestral) class through the polymorphic bins to the
# fixed-derived (divergence) class in the last entry
sfs_neut = fd.Spectrum([171150, 997, 441, 228, 156, 117, 114, 83, 105, 109, 6500])
sfs_sel = fd.Spectrum([793221, 1329, 499, 265, 162, 104, 117, 90, 94, 119, 14000])

# specifying the divergence target sizes includes divergence in the likelihood
inf = fd.BaseInference(
    sfs_neut=sfs_neut,
    sfs_sel=sfs_sel,
    n_sites_div_neut=180000,
    n_sites_div_sel=810000,
    fixed_params=dict(all=dict(eps=0, h=0.5))
)

inf.run()

inf.plot_discretized();

INFO:BaseInference: Using divergence counts in the likelihood.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [--:--<--:--, --it/s]
INFO:Optimization: Optimizing 4 parameters: [all.S_d, all.b, all.p_b, all.S_b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.p_b: (0, 0.0044950944, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1.059e+04 ± 2.8e+04, all.b: 0.1545 ± 3.1, all.p_b: 0.004495 ± 0.13, all.S_b: 63.22 ± 24, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -45.75 ± 1.6e+02} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.293e+04 ± 3.3e+04, all.b: 0.1698 ± 0.039, all.p_b: 0.008446 ± 0.0088, all.S_b: 65.34 ± 36, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -52 ± 5, i_best_run: 0.51 ± 0.5, likelihoods_std: 7.625 ± 34} (mean ± std)

../_images/5ae1f835a56957c676c3b88c27a0286b08bac09dd7b8f1a846b7dbf8ff9136b1.png

R

# the SFS runs from the monomorphic (ancestral) class through the polymorphic bins to the
# fixed-derived (divergence) class in the last entry
sfs_neut <- fd$Spectrum(c(171150, 997, 441, 228, 156, 117, 114, 83, 105, 109, 6500))
sfs_sel <- fd$Spectrum(c(793221, 1329, 499, 265, 162, 104, 117, 90, 94, 119, 14000))

# specifying the divergence target sizes includes divergence in the likelihood
inf <- fd$BaseInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel,
  n_sites_div_neut = 180000,
  n_sites_div_sel = 810000,
  fixed_params = list(all = list(eps = 0, h = 0.5))
)

sfs_modelled <- inf$run()

p <- inf$plot_discretized()

INFO:BaseInference: Using divergence counts in the likelihood.
INFO:Discretization: Precomputing semidominant DFE-SFS transformation using midpoint integration.
Discretization>Precomputing: 100%|██████████| 9/9 [--:--<--:--, --it/s]
INFO:Optimization: Optimizing 4 parameters: [all.S_d, all.b, all.p_b, all.S_b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]
WARNING:Optimization: The MLE estimate is close to the upper bound for {} and lower bound for {all.p_b: (0, 0.0044950944, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
INFO:BaseInference: Inference results: {all.S_d: -1.059e+04 ± 2.8e+04, all.b: 0.1545 ± 3.1, all.p_b: 0.004495 ± 0.13, all.S_b: 63.22 ± 24, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -45.75 ± 1.6e+02} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.293e+04 ± 3.3e+04, all.b: 0.1698 ± 0.039, all.p_b: 0.008446 ± 0.0088, all.S_b: 65.34 ± 36, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -52 ± 5, i_best_run: 0.51 ± 0.5, likelihoods_std: 7.625 ± 34} (mean ± std)

../_images/8ca8db18fa50b981f9acad3f5afdccc8a0cb8f00b6d895b44e536c21a37411b1.png

The estimator used by get_alpha() follows the inference mode, but can be switched explicitly via its use_divergence argument, which allows \(\alpha\) estimated from divergence to be compared with \(\alpha\) estimated from polymorphism alone. get_omega() returns \(\omega\), the ratio of non-synonymous to synonymous substitution rates (\(d_N/d_S\)), and get_omega_a() returns its adaptive component \(\omega_a = \alpha\,\omega\).

Python

print(f'alpha (with divergence):   {inf.get_alpha():.7g}')
print(f'alpha (polymorphism only): {inf.get_alpha(use_divergence=False):.7g}')
print(f'omega:   {inf.get_omega():.7g}')
print(f'omega_a: {inf.get_omega_a():.7g}')

alpha (with divergence):   0.6614585
alpha (polymorphism only): 0.5940308
omega:   0.4786325
omega_a: 0.3165955

R

cat('alpha (with divergence):  ', inf$get_alpha(), '\n')
cat('alpha (polymorphism only):', inf$get_alpha(use_divergence = FALSE), '\n')
cat('omega:  ', inf$get_omega(), '\n')
cat('omega_a:', inf$get_omega_a(), '\n')

alpha (with divergence):   0.6614585
alpha (polymorphism only): 0.5940308
omega:   0.4786325
omega_a: 0.3165955

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, whereas in practice errors may differ across classes. Nevertheless, below, we infer the full DFE while allowing eps to vary.

Python

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

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 5 parameters: [all.S_d, all.b, all.p_b, all.S_b, all.eps].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.064e+04 ± 3.6e+03, all.b: 0.1508 ± 0.0036, all.p_b: 0 ± 0.0018, all.S_b: 0.0001247 ± 0.051, all.eps: 0.006854 ± 0.00084, all.h: 0.5 ± 0, likelihood: -34.63 ± 0.052} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.256e+04 ± 3.5e+04, all.b: 0.9817 ± 2.5, all.p_b: 0.06955 ± 0.093, all.S_b: 0.01596 ± 0.093, all.eps: 0.006245 ± 0.0071, all.h: 0.5 ± 0, likelihood: -40.12 ± 4.6, i_best_run: 0.54 ± 0.5, likelihoods_std: 0.04447 ± 0.3} (mean ± std)

R

inf <- fd$BaseInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel,
  fixed_params = list(all = list(h = 0.5))
)

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 5 parameters: [all.S_d, all.b, all.p_b, all.S_b, all.eps].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.064e+04 ± 3.6e+03, all.b: 0.1508 ± 0.0036, all.p_b: 0 ± 0.0018, all.S_b: 0.0001247 ± 0.051, all.eps: 0.006854 ± 0.00084, all.h: 0.5 ± 0, likelihood: -34.63 ± 0.052} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -2.256e+04 ± 3.5e+04, all.b: 0.9817 ± 2.5, all.p_b: 0.06955 ± 0.093, all.S_b: 0.01596 ± 0.093, all.eps: 0.006245 ± 0.0071, all.h: 0.5 ± 0, likelihood: -40.12 ± 4.6, i_best_run: 0.54 ± 0.5, likelihoods_std: 0.04447 ± 0.3} (mean ± std)

The bootstrap distribution of eps:

Python

inf.bootstraps.eps.hist(grid=False, figsize=(4.8, 2.5));

../_images/ce0ca8e13c54229c4c1942ee5fe39242f970ddcf931e55c562c5dc3988ac0e3f.png

R

par(mar = c(2.5, 3, 1, 1))
hist(inf$bootstraps$eps, main = "", xlab = "", col = "#1f77b4", border = "white")

../_images/55fc5b991a3fc719c7da19612b775e7bf1d725e0d0cd7de3d4e0d28d99a3764c.png

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 does not significantly improve the model fit. We can check this in a more principled way by performing a likelihood-ratio test as done below.

Nested model comparison#

The significance of including ancestral-allele misidentification and beneficial mutations can be assessed with likelihood ratio tests, using plot_nested_models(). The LRTs compare the likelihood of the inferred DFE to the likelihood of a nested model where some parameters are held fixed. Alternatively, compare_nested() directly compares two nested models.

Python

# 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 [--:--<--:--, --it/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.0043298919, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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 [--:--<--:--, --it/s]
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]

../_images/b7ced19e507132f153c30aa243e1ea0d9f50ff0a95565aab8b9dd51335efcb07.png

R

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

p <- inf$plot_nested_models()

fd$logger$setLevel('INFO')

BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.0043298919, 0.5)} [(lower, value, upper)], but this might be nothing to worry about.
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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 [--:--<--:--, --it/s]
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]

../_images/e410a31e5280da50e2c2bfaa30e5d39b793a053c3cf8702bf4262dd43a770fde.png

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

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 practice, 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 scaled selection coefficient S = 4 Ne s.

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.

Python

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:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:Discretization: Precomputing DFE-SFS transformation for fixed dominance coefficients.
Discretization>Precomputing: 100%|██████████| 1809/1809 [--:--<--:--, --it/s]
INFO:Optimization: Optimizing 2 parameters: [all.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.159. 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.1411 ± 2.6e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0 ± 0, likelihood: -75.43 ± 1.2e-12} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -1e+05 ± 0, all.b: 0.1404 ± 0.0035, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0 ± 0, likelihood: -82.36 ± 15, i_best_run: 0.53 ± 0.5, likelihoods_std: 0.06262 ± 0.53} (mean ± std)

R

inf <- fd$BaseInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel,
  fixed_params = list(all = list(eps = 0, h = 0, p_b = 0, S_b = 1)),
  h_callback = function(h, S) 0.4 * exp(-0.1 * abs(S))
)

sfs_modelled <- inf$run()

INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:Discretization: Precomputing DFE-SFS transformation for fixed dominance coefficients.
Discretization>Precomputing: 100%|██████████| 1809/1809 [--:--<--:--, --it/s]
INFO:Optimization: Optimizing 2 parameters: [all.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.159. 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.1411 ± 2.6e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0 ± 0, likelihood: -75.43 ± 1.2e-12} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -1e+05 ± 0, all.b: 0.1404 ± 0.0035, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0 ± 0, likelihood: -82.36 ± 15, i_best_run: 0.53 ± 0.5, likelihoods_std: 0.06262 ± 0.53} (mean ± std)

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

Python

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

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

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.621e+04 ± 3.8, all.b: 0.1305 ± 1.5e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 4.6e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.094e+04 ± 4.1e+04, all.b: 0.137 ± 0.022, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.59 ± 5.3, i_best_run: 0.49 ± 0.5, likelihoods_std: 0.06163 ± 0.57} (mean ± std)

../_images/52a7b7e79bcfe24d3fba3bb8422106a7b8c30fc083dfdb8176f52bc542cb7d85.png

R

inf2 <- fd$BaseInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel,
  fixed_params = list(all = list(eps = 0, h = 0.5, p_b = 0, S_b = 1))
)

sfs_modelled <- inf2$run()

p <- fd$DFE$plot_many(list(inf$get_dfe(), inf2$get_dfe()), labels = c('partly recessive', 'h=0.5'))

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.621e+04 ± 3.8, all.b: 0.1305 ± 1.5e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 4.6e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -5.094e+04 ± 4.1e+04, all.b: 0.137 ± 0.022, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -42.59 ± 5.3, i_best_run: 0.49 ± 0.5, likelihoods_std: 0.06163 ± 0.57} (mean ± std)

../_images/85efdfd533021b716fdb0c9d8aff6f6cd32c297bc6ca7b3f22c2fc7d3ae3b636.png

Assuming that 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, folded spectra are passed 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.

Python

import matplotlib.pyplot as plt

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

inf.run()

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

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

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: -1.654e+04 ± 0.87, all.b: 0.1464 ± 9.4e-07, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -21.58 ± 1.3e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -3.885e+04 ± 4e+04, all.b: 0.1518 ± 0.032, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -24.88 ± 3.5, i_best_run: 0.49 ± 0.5, likelihoods_std: 0.003806 ± 0.022} (mean ± std)

../_images/f78a8753a73605c7b515f13e904779c5a5d9995f4ec334f748ec0983d7cdd8e0.png

R

inf <- fd$BaseInference(
  sfs_neut = sfs_neut$fold(),
  sfs_sel = sfs_sel$fold()
)

sfs_modelled <- inf$run()

# plot the inferred DFE and the SFS comparison
p1 <- inf$plot_discretized(show = FALSE, intervals = c(-Inf, -100, -10, -1, 0))
p2 <- inf$plot_sfs_comparison(show = FALSE)

cowplot::plot_grid(p1, p2, ncol = 2)

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: -1.654e+04 ± 0.87, all.b: 0.1464 ± 9.4e-07, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -21.58 ± 1.3e-09} (best_run ± std_across_runs)
BaseInference>Bootstrapping (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -3.885e+04 ± 4e+04, all.b: 0.1518 ± 0.032, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -24.88 ± 3.5, i_best_run: 0.49 ± 0.5, likelihoods_std: 0.003806 ± 0.022} (mean ± std)

../_images/70b0b2af6bbc94182e706ba456942042b752c611b94e57c27519a1c5746735ca.png

Serialization#

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

Python

# save the inference object to a file, which BaseInference.from_file restores
inf.to_file("serialized.json")

# save a short summary to a file
inf.get_summary().to_file("summary.json")

R

# save the inference object to a file, which BaseInference$from_file restores
inf$to_file("serialized.json")

# save a short summary to a file
inf$get_summary()$to_file("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 that share S_d, the mean selection coefficient of deleterious mutations (cf. GammaExpParametrization). For more complex stratifications, see the Parser module.

Python

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

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

inf.run();

INFO:JointInference: No divergence counts provided, inferring from polymorphism only.
INFO:JointInference: Using shared parameters [SharedParams(params=['S_d'], types=['pendula', 'pubescens'])].
INFO:JointInference: Including covariates: {}.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:JointInference: Running marginal inference for type 'all'.
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]
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.1e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -45.09 ± 5.3e-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.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 3.2, all.b: 0.1305 ± 1.3e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 2.8e-09} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'pubescens'.
INFO:Optimization: Optimizing 2 parameters: [all.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.1035 ± 1.6e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -43.81 ± 4.5e-12} (best_run ± std_across_runs)
INFO:JointInference: Running joint inference for types ['pendula', 'pubescens'].
INFO:Optimization: Optimizing 3 parameters: [pendula.b, pubescens.b, pendula:pubescens.S_d].
JointInference>Performing joint inference: 100%|██████████| 10/10 [--:--<--:--, --it/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 ± 5.6e-09, pendula.p_b: 0 ± 0, pendula.S_b: 1 ± 0, pendula.eps: 0 ± 0, pendula.h: 0.5 ± 0, pubescens.b: 0.1035 ± 7.5e-09, 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 ± 0, likelihood: -79.46 ± 7.7e-11} (best_run ± std_across_runs)
INFO:JointInference: Bootstrapping type 'all'.
BaseInference>Bootstrapping 'all' (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -7.935e+04 ± 3e+04, all.b: 0.1105 ± 0.0073, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -53.76 ± 5.9, i_best_run: 0.53 ± 0.5, likelihoods_std: 3.636 ± 21} (mean ± std)
INFO:JointInference: Bootstrapping type 'pendula'.
BaseInference>Bootstrapping 'pendula' (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)
INFO:JointInference: Bootstrapping type 'pubescens'.
BaseInference>Bootstrapping 'pubescens' (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -9.105e+04 ± 2.1e+04, all.b: 0.1048 ± 0.005, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -52.64 ± 6.7, i_best_run: 0.4 ± 0.49, likelihoods_std: 0.6082 ± 6} (mean ± std)
JointInference>Bootstrapping joint inference (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:JointInference: Bootstrap summary: {pendula.b: 0.1213 ± 0.0088, pendula.p_b: 0 ± 0, pendula.S_b: 1 ± 0, pendula.eps: 0 ± 0, pendula.h: 0.5 ± 0, pendula.S_d: -8.33e+04 ± 2.9e+04, pubescens.b: 0.1068 ± 0.007, pubescens.p_b: 0 ± 0, pubescens.S_b: 1 ± 0, pubescens.eps: 0 ± 0, pubescens.h: 0.5 ± 0, pubescens.S_d: -8.33e+04 ± 2.9e+04, likelihood: -96.19 ± 9.6, i_best_run: 0.59 ± 0.49, likelihoods_std: 0.01092 ± 0.078} (mean ± std)

R

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

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

inf <- fd$JointInference(
  sfs_neut = sfs_neut,
  sfs_sel = sfs_sel,
  shared_params = list(fd$SharedParams(types = c("pendula", "pubescens"), params = list("S_d")))
)

sfs_modelled <- inf$run()

INFO:JointInference: No divergence counts provided, inferring from polymorphism only.
INFO:JointInference: Using shared parameters [SharedParams(params=['S_d'], types=['pendula', 'pubescens'])].
INFO:JointInference: Including covariates: {}.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:BaseInference: No divergence counts provided, inferring from polymorphism only.
INFO:JointInference: Running marginal inference for type 'all'.
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]
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.1e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -45.09 ± 5.3e-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.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/s]
INFO:BaseInference: Inference results: {all.S_d: -3.389e+04 ± 3.2, all.b: 0.1305 ± 1.3e-06, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -35.44 ± 2.8e-09} (best_run ± std_across_runs)
INFO:JointInference: Running marginal inference for type 'pubescens'.
INFO:Optimization: Optimizing 2 parameters: [all.S_d, all.b].
BaseInference>Performing inference: 100%|██████████| 10/10 [--:--<--:--, --it/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.1035 ± 1.6e-09, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -43.81 ± 4.5e-12} (best_run ± std_across_runs)
INFO:JointInference: Running joint inference for types ['pendula', 'pubescens'].
INFO:Optimization: Optimizing 3 parameters: [pendula.b, pubescens.b, pendula:pubescens.S_d].
JointInference>Performing joint inference: 100%|██████████| 10/10 [--:--<--:--, --it/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 ± 5.6e-09, pendula.p_b: 0 ± 0, pendula.S_b: 1 ± 0, pendula.eps: 0 ± 0, pendula.h: 0.5 ± 0, pubescens.b: 0.1035 ± 7.5e-09, 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 ± 0, likelihood: -79.46 ± 7.7e-11} (best_run ± std_across_runs)
INFO:JointInference: Bootstrapping type 'all'.
BaseInference>Bootstrapping 'all' (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -7.935e+04 ± 3e+04, all.b: 0.1105 ± 0.0073, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -53.76 ± 5.9, i_best_run: 0.53 ± 0.5, likelihoods_std: 3.636 ± 21} (mean ± std)
INFO:JointInference: Bootstrapping type 'pendula'.
BaseInference>Bootstrapping 'pendula' (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)
INFO:JointInference: Bootstrapping type 'pubescens'.
BaseInference>Bootstrapping 'pubescens' (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:BaseInference: Bootstrap summary: {all.S_d: -9.105e+04 ± 2.1e+04, all.b: 0.1048 ± 0.005, all.p_b: 0 ± 0, all.S_b: 1 ± 0, all.eps: 0 ± 0, all.h: 0.5 ± 0, likelihood: -52.64 ± 6.7, i_best_run: 0.4 ± 0.49, likelihoods_std: 0.6082 ± 6} (mean ± std)
JointInference>Bootstrapping joint inference (2 runs each): 100%|██████████| 100/100 [--:--<--:--, --it/s]
INFO:JointInference: Bootstrap summary: {pendula.b: 0.1213 ± 0.0088, pendula.p_b: 0 ± 0, pendula.S_b: 1 ± 0, pendula.eps: 0 ± 0, pendula.h: 0.5 ± 0, pendula.S_d: -8.33e+04 ± 2.9e+04, pubescens.b: 0.1068 ± 0.007, pubescens.p_b: 0 ± 0, pubescens.S_b: 1 ± 0, pubescens.eps: 0 ± 0, pubescens.h: 0.5 ± 0, pubescens.S_d: -8.33e+04 ± 2.9e+04, likelihood: -96.19 ± 9.6, i_best_run: 0.59 ± 0.49, likelihoods_std: 0.01092 ± 0.078} (mean ± std)

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

Python

inf.plot_inferred_parameters();

../_images/2f8fc94384cb49811ace52fe3da53fccc76a513afe0750a139508a7d70ac2aa1.png

R

p <- inf$plot_inferred_parameters()

../_images/fa351e958a4665d5fd6eb16ae3315b1e489b68715af525ae07971aaec68542ad.png

marginal.pendula and marginal.pubescens are the marginal inferences for the respective type. marginal.all is the marginal inference obtained by adding up the spectra of all types. joint.pendula and joint.pubescens are the joint inferences for the respective type. We can see that S_d is 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 is available through marginal_inferences.

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

Python

inf.plot_discretized();

../_images/b3f7671677ae5d4d736b13bd1ff5181d68787fbbbb287a461f476d4868aa10ab.png

R

p <- inf$plot_discretized()

../_images/806cb8ad353e5e87c2cbc33b03f3b49ee91346e8e09baff8f35a58770d4e1175.png

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.

Python

inf.perform_lrt_shared()

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

0.5097493

R

inf$perform_lrt_shared()

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

[1] 0.5097493

The test is not significant, indicating that the simpler model of sharing the parameters explains the data sufficiently well. Indeed, the inferred parameters of the joint and the marginal inferences differ little.