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

stdpopsim/catalog/DroMel/dfes.py

@@ -47,14 +47,74 @@ def _HuberDFE():
         description=description,
         long_description=long_description,
         mutation_types=[neutral, negative],
-        proportions=[prop_synonymous, 1 - prop_synonymous],  # LNS = 2.85 x LS
+        proportions=[prop_synonymous, 1 - prop_synonymous],
         citations=citations,
     )
 
 
 _species.add_dfe(_HuberDFE())
 
 
+def _ZhenDFE():
+    id = "GammaPos_H17"
+    description = "Deleterious Gamma DFE with fixed-s beneficials"
+    long_description = """
+    The DFE estimated from D.melanogaster-simulans data estimated in Zhen et al
+    (2021, https://dx.doi.org/10.1101/gr.256636.119).
+    This uses the demographic model and deleterious-only DFE from
+    Huber et al (2017), and then the number of nonsynonymous differences
+    to simulans to estimate a proportion of nonsynonmyous differences
+    and (single) selection coefficient. So, this is the "Gamma_H17" DFE,
+    with some proportion of positive selection with fixed s.
+    """
+    citations = [
+        stdpopsim.Citation(
+            author="Zhen et al.",
+            year=2021,
+            doi="https://dx.doi.org/10.1101/gr.256636.119",
+            reasons={stdpopsim.CiteReason.DFE},
+        )
+    ]
+    neutral = stdpopsim.MutationType()
+    gamma_shape = 0.33  # shape
+    gamma_mean = -3.96e-04  # expected value
+    h = 0.5  # dominance coefficient
+    negative = stdpopsim.MutationType(
+        dominance_coeff=h,
+        distribution_type="g",  # gamma distribution
+        distribution_args=[gamma_mean, gamma_shape],
+    )
+    # p. 2 in supplement says that the total sequence length of synonymous sites LS
+    # related to the total sequence length of nonsynonymous sites LNS
+    # by LNS = 2.85 * LS
+    # so this is 1 / (1 + 2.85) = 0.2597402597402597
+    prop_synonymous = 0.26
+
+    sel_coeff = 10 ** (-4.801)
+    prop_beneficials = 6.75e-4
+    positive = stdpopsim.MutationType(
+        dominance_coeff=0.5,
+        distribution_type="f",
+        distribution_args=[sel_coeff],
+    )
+
+    return stdpopsim.DFE(
+        id=id,
+        description=description,
+        long_description=long_description,
+        mutation_types=[neutral, negative, positive],
+        proportions=[
+            prop_synonymous,
+            (1 - prop_synonymous) * (1 - prop_beneficials),
+            (1 - prop_synonymous) * prop_beneficials,
+        ],
+        citations=citations,
+    )
+
+
+_species.add_dfe(_ZhenDFE())
+
+
 def _RagsdaleDFE():
     id = "LognormalPlusPositive_R16"
     description = "Deleterious log-normal and beneficial mixed DFE"

stdpopsim/catalog/HomSap/dfes.py

@@ -259,3 +259,66 @@ def _RodriguesDFE():
 
 
 _species.add_dfe(_RodriguesDFE())
+
+
+def _ZhenDFE():
+    id = "GammaPos_Z21"
+    description = "Deleterious Gamma DFE with fixed-s beneficials"
+    long_description = """
+    The DFE estimated from human-chimp data estimated in Zhen et al
+    (2021, https://dx.doi.org/10.1101/gr.256636.119).
+    This uses the demographic model and deleterious-only DFE from
+    Huber et al (2017), and then the number of nonsynonymous differences
+    to chimpanzee to estimate a proportion of nonsynonmyous differences
+    and (single) selection coefficient. So, this is the "Gamma_H17" DFE,
+    with some proportion of positive selection with fixed s.
+    """
+    citations = [
+        stdpopsim.Citation(
+            author="Zhen et al.",
+            year=2021,
+            doi="https://dx.doi.org/10.1101/gr.256636.119",
+            reasons={stdpopsim.CiteReason.DFE},
+        )
+    ]
+    # modified from _HuberDFE() above
+    neutral = stdpopsim.MutationType()
+    gamma_shape = 0.19  # shape
+    # Extra factor of 2 in mean is to account for difference between tools
+    # in how fitness is defined
+    # (1+s for homozygote in SLiM versus 1+2s in dadi)
+    gamma_mean = -0.014 * 2  # expected value
+    negative = stdpopsim.MutationType(
+        dominance_coeff=0.5,
+        distribution_type="g",  # gamma distribution
+        distribution_args=[gamma_mean, gamma_shape],
+    )
+    # p. 2 in supplement says that the total sequence length of synonymous sites LS
+    # related to the total sequence length of nonsynonymous sites LNS
+    # by LNS = 2.31 * LS
+    # so, this is 1 / (1 + 2.31) = 0.3021148036253776
+    prop_synonymous = 0.3
+
+    sel_coeff = 10 ** (-3.949)
+    prop_beneficial = 1.55e-2
+    positive = stdpopsim.MutationType(
+        dominance_coeff=0.5,
+        distribution_type="f",
+        distribution_args=[sel_coeff],
+    )
+
+    return stdpopsim.DFE(
+        id=id,
+        description=description,
+        long_description=long_description,
+        mutation_types=[neutral, negative, positive],
+        proportions=[
+            prop_synonymous,
+            (1 - prop_synonymous) * (1 - prop_beneficial),
+            (1 - prop_synonymous) * prop_beneficial,
+        ],
+        citations=citations,
+    )
+
+
+_species.add_dfe(_ZhenDFE())