diff --git a/stdpopsim/catalog/AraTha/dfes.py b/stdpopsim/catalog/AraTha/dfes.py index fbcbdd092..fad3117ad 100644 --- a/stdpopsim/catalog/AraTha/dfes.py +++ b/stdpopsim/catalog/AraTha/dfes.py @@ -10,19 +10,15 @@ def _HuberDFE(): - id = "Gamma_H18" - description = "Deleterious Gamma DFE" + id = "GammaAdditive_H18" + description = "Deleterious additive Gamma DFE" long_description = """ - Return negative (no neutral, in this case) MutationType()s - representing an Arabidopsis DFE. From Huber et al. (2018). - Gamma parameters are based on Supplementary Table 4, - the genome-wide, additive-only model for A. LYRATA- due to - challenges with simulating with selfing for A. thaliana. - The Supplementary Table 4 DFEs are not noted to contain - any neutral mutations (in contrast, Supplementary Table 3 - notes neutral proportions for two other DFEs), and in the - main text including neutral mutations seems like a - supplementary analysis rather than the main strategy. + Additive, deleterious DFE from Huber et al. (2018) for Arabidopsis, + estimated from the SFS of A. lyrata as a Gamma distribution + of deleterious effects. Parameters are from Supplementary Table 4, + the "genome-wide, additive-only model for A. lyrata". The DFE for + A. lyrata (rather than A. thaliana) is provided due to challenges + with simulating selfing. """ citations = [ stdpopsim.Citation( @@ -34,12 +30,15 @@ def _HuberDFE(): ] neutral = stdpopsim.MutationType() gamma_shape = 0.27 # shape - gamma_mean = -0.0004 # expected value + gamma_scale = -0.0004 # scale + gamma_mean = gamma_shape * gamma_scale h = 0.5 # dominance coefficient negative = stdpopsim.MutationType( dominance_coeff=h, distribution_type="g", # gamma distribution - distribution_args=[gamma_mean, gamma_shape], + # extra factor of 2 is to convert dadi to SLiM + # (1+s for homozygote in SLiM versus 1+2s in dadi) + distribution_args=[-2 * gamma_mean, gamma_shape], ) return stdpopsim.DFE( diff --git a/stdpopsim/qc/AraTha.py b/stdpopsim/qc/AraTha.py index 6fd5ec96f..8ef64d3a5 100644 --- a/stdpopsim/qc/AraTha.py +++ b/stdpopsim/qc/AraTha.py @@ -208,3 +208,29 @@ def HuberThreeEpoch(): _species.get_demographic_model("African3Epoch_1H18").register_qc(HuberThreeEpoch()) + + +def HuberDFE(): + id = "QC-GammaAdditive_H18" + # Huber et al. 2018 Supplementary table 4. + # genome-wide A. lyrata additive fitness effects + neutral = stdpopsim.MutationType() + gamma_shape = 0.27 + gamma_scale = 0.0004 + gamma_mean = -(gamma_shape * gamma_scale) # negative since deleterious mutations + h = 0.5 + negative = stdpopsim.MutationType( + dominance_coeff=h, + distribution_type="g", + distribution_args=[-2 * gamma_mean, gamma_shape], + ) + return stdpopsim.DFE( + id=id, + description=id, + long_description=id, + mutation_types=[neutral, negative], + proportions=[0.0, 1.0], + ) + + +_species.get_dfe("GammaAdditive_H18").register_qc(HuberDFE()) diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py index 722dec002..810285203 100644 --- a/tests/test_slim_engine.py +++ b/tests/test_slim_engine.py @@ -1213,7 +1213,7 @@ def test_not_simulated_outside_region(self): left, right = 100000, 900000 contig = species.get_contig("1", left=left, right=right) - dfe = species.get_dfe("Gamma_H18") + dfe = species.get_dfe("GammaAdditive_H18") exons = species.get_annotations("araport_11_exons") exon_intervals = exons.get_chromosome_annotations("1") contig.add_dfe(intervals=exon_intervals, DFE=dfe)