From cb1550ab1a59f90e5034c597eceab9de4555945d Mon Sep 17 00:00:00 2001
From: Clara Rehmann <66445895+clararehmann@users.noreply.github.com>
Date: Wed, 22 Jan 2025 20:56:00 -0800
Subject: [PATCH] Dromel qc (#1668)

* qc for DroMel DFE
---
 stdpopsim/catalog/DroMel/dfes.py | 13 +++++++++----
 stdpopsim/qc/DroMel.py           | 30 ++++++++++++++++++++++++++++++
 2 files changed, 39 insertions(+), 4 deletions(-)

diff --git a/stdpopsim/catalog/DroMel/dfes.py b/stdpopsim/catalog/DroMel/dfes.py
index d7c8cada4..0d3222ef1 100644
--- a/stdpopsim/catalog/DroMel/dfes.py
+++ b/stdpopsim/catalog/DroMel/dfes.py
@@ -29,20 +29,25 @@ def _HuberDFE():
     ]
     neutral = stdpopsim.MutationType()
     gamma_shape = 0.33  # shape
-    gamma_mean = -3.96e-04  # expected value
+    gamma_scale = 1.2e-3  # scale
+    gamma_mean = gamma_shape * gamma_scale  # expected value
     h = 0.5  # dominance coefficient
     negative = stdpopsim.MutationType(
         dominance_coeff=h,
         distribution_type="g",  # gamma distribution
-        distribution_args=[gamma_mean, gamma_shape],
+        # (1+s for homozygote in SLiM versus 1+2s in dadi)
+        distribution_args=[-2 * 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
     return stdpopsim.DFE(
         id=id,
         description=description,
         long_description=long_description,
         mutation_types=[neutral, negative],
-        proportions=[0.26, 0.74],  # LNS = 2.85 x LS
+        proportions=[prop_synonymous, 1 - prop_synonymous],  # LNS = 2.85 x LS
         citations=citations,
     )
 
diff --git a/stdpopsim/qc/DroMel.py b/stdpopsim/qc/DroMel.py
index f7458b25a..843033513 100644
--- a/stdpopsim/qc/DroMel.py
+++ b/stdpopsim/qc/DroMel.py
@@ -115,3 +115,33 @@ def SheehanSongThreeEpic():
 
 
 _species.get_demographic_model("African3Epoch_1S16").register_qc(SheehanSongThreeEpic())
+
+
+def Gamma_H17():
+    id = "Gamma_H17"
+    description = "Deleterious Gamma DFE"
+    long_description = "QC version"
+    neutral = stdpopsim.MutationType()
+    gamma_shape = 0.33
+    gamma_scale = 1.2e-3
+    gamma_mean = gamma_shape * gamma_scale
+    h = 0.5  # dominance coefficient
+    negative = stdpopsim.MutationType(
+        dominance_coeff=h,
+        distribution_type="g",  # gamma distribution
+        # (1+s for homozygote in SLiM versus 1+2s in dadi)
+        distribution_args=[-2 * gamma_mean, gamma_shape],
+    )
+    # LNS = 2.85 * LS
+    # prop_synonymous = 1/(1+2.85) = 0.26
+    prop_synonymous = 0.26
+    return stdpopsim.DFE(
+        id=id,
+        description=description,
+        long_description=long_description,
+        mutation_types=[neutral, negative],
+        proportions=[prop_synonymous, 1 - prop_synonymous],
+    )
+
+
+_species.get_dfe("Gamma_H17").register_qc(Gamma_H17())