almost finished

stdpopsim/catalog/HomSap/dfes.py

@@ -161,42 +161,46 @@ def _KyriazisDFE():
     0.63 is assumed, informed by analysis of 30Mb of human coding sequences.
     The model is also augmented with an additional proportion of 0.3% of
     recessive lethals, based on the analysis of Wade et al. (2023).
-#    Gamma parameters are given as the shape and the mean selection coefficient
- #   (E[s]) escaled, such that E[s] = -shape*scale*2/(2Na),
-  #  # where Na is the ancestral population size (12378) as in Kim et al. (2017).
     """
     citations = [
         stdpopsim.Citation(
             author="Kyriazis et al.",
             year=2023,
             doi="https://doi.org/10.1086/726736",
-#            reasons={stdpopsim.CiteReason.DFE}, ?
+            reasons={stdpopsim.CiteReason.DFE},
         )
     ]
     neutral = stdpopsim.MutationType()
     gamma_mean = 0.0131
     gamma_shape = 0.186
-    # 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
     coeffs = [0.45,  0.2,  0.05, 0]
     breaks = [0.001, 0.01, 0.1, 1]
-    negative = stdpopsim.MutationType(
+    gamma = stdpopsim.MutationType(
         dominance_coeff_list = coefs,
         dominance_coeff_breaks = breaks,
         distribution_type="g",  # gamma distribution
         distribution_args=[gamma_mean, gamma_shape],
     )
-
+    lethal = stdpopsim.MutationType(
+        dominance_coeff_list = coefs,
+        dominance_coeff_breaks = breaks,
+        <<selection_coeff???>> = 1      ## *******
+    )    
+    # To get 0.63 deleterious mutations per 30Mb diploid coding sequence with 
+    # mu=1.29e-8 (current mu from catalog), we need a deleterious proportion
+    # of 0.8139535. Next, 0.03% of deleterious are lethals:
+    gamma_prop = 0.8139535 * 0.997
+    lethal_prop = 0.8139535 * 0.003
+    neutral_prop = 1 - gamma_prop - lethal_prop
+
     # To get 0.63 deleterious mutations per 30Mb diploid coding sequence with
     # mu=1.29e-8 (current mu from catalog), we need a proportion of 0.8139535:
     return stdpopsim.DFE(
         id=id,
         description=description,
         long_description=long_description,
-        mutation_types=[neutral, negative],
-        proportions=[0.1860465, 0.8139535],
+        mutation_types=[neutral, gamma, lethal],
+        proportions=[neutral_prop, gamma_prop, lethal_prop],
         citations=citations,
     )