positive-containing DFE from Zhen et al 2021 (GammaPos_Z21, HomSap and DroMel)

[petrelharp]

This is as recommended by @klohmueller over in #1469. It's really just Huber et al's Gamma_H17 except a proportion p of the changes are beneficial with coefficient s; the paper shows that p and s are highly confounded (and ps can be estimated better); but these are the MLE values. These values are using divergence to chimp.

Kirk said:

For a DFE with positive selection in humans—we could use this one from Table S6 of Zhen et al. (2021). This DFE is for nonsynonymous mutations. We inferred 1.55% of nonsynonymous mutations were positively selected in humans with a s of 0.0001124605. The remaining nonsynonymous mutations came from a DFE for negative selection (parameters for this are in Table S3). The shape parameter is 0.19 and the scale (beta) in terms of s is 0.074.

What I've implemented is this, and I've verified all those parameters above - so, I figure this is actually mostly QC'ed (since Kirk pulled out those numbers and I verified them). So, I'm tempted to just put in the QC code also - but that wouldn't catch coding errors.

Edit: I've also added the estimated one for DroMelfrom Table S6 as above.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.85%. Comparing base (038c506) to head (4780465).
Report is 2 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #1656   +/-   ##
=======================================
  Coverage   99.85%   99.85%           
=======================================
  Files         139      139           
  Lines        4743     4774   +31     
  Branches      470      470           
=======================================
+ Hits         4736     4767   +31     
  Misses          3        3           
  Partials        4        4           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[silastittes]

ran fine after adding prop_synonymous = 0.26 back to line 44.

Here's what I ran locally:

import stdpopsim
import numpy as np

species = stdpopsim.get_species("DroMel")
dfe = species.get_dfe("GammaPos_H17")
contig = species.get_contig(length=3000)
model = species.get_demographic_model("OutOfAfrica_2L06")
samples = {"AFR": 50, "EUR": 50}

gene_interval = np.array([[10000, 20000]])
contig.add_dfe(intervals=gene_interval, DFE=dfe)

engine = stdpopsim.get_engine("slim")
ts = engine.simulate(
    model,
    contig,
    samples,
    seed=236,
    slim_scaling_factor=1000,
    slim_burn_in=1,
)

print(ts)

║TreeSequence               ║
╠═══════════════╤═══════════╣
║Trees          │       2200║
╟───────────────┼───────────╢
║Sequence Length│       3000║
╟───────────────┼───────────╢
║Time Units     │generations║
╟───────────────┼───────────╢
║Sample Nodes   │        200║
╟───────────────┼───────────╢
║Total Size     │  793.0 KiB║
╚═══════════════╧═══════════╝
╔═══════════╤══════╤═════════╤════════════╗
║Table      │Rows  │Size     │Has Metadata║
╠═══════════╪══════╪═════════╪════════════╣
║Edges      │14,380│449.4 KiB│          No║
╟───────────┼──────┼─────────┼────────────╢
║Individuals│   100│ 11.6 KiB│         Yes║
╟───────────┼──────┼─────────┼────────────╢
║Migrations │     0│  8 Bytes│          No║
╟───────────┼──────┼─────────┼────────────╢
║Mutations  │ 1,065│ 57.6 KiB│         Yes║
╟───────────┼──────┼─────────┼────────────╢
║Nodes      │ 2,733│ 97.6 KiB│         Yes║
...
╟───────────┼──────┼─────────┼────────────╢
║Sites      │   894│ 21.8 KiB│          No║
╚═══════════╧══════╧═════════╧════════════╝

[silastittes]

Updates worked for me!