Dromel qc

[clararehmann]

Addresses #1054

I used the gamma distribution parameters from the full model (no singletons, recent mutation rate estimates) in table S2 and calculated the proportions of neutral and gamma from Peter's comment:

Note that the proportions of neutral & gamma ought to come from this snippet (see this explanation in #1100).

(I have not checked.)

Originally posted by @petrelharp in #1054

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.85%. Comparing base (5e63e81) to head (258bdf1).
Report is 1 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #1668   +/-   ##
=======================================
  Coverage   99.84%   99.85%           
=======================================
  Files         136      136           
  Lines        4656     4671   +15     
  Branches      470      470           
=======================================
+ Hits         4649     4664   +15     
  Misses          3        3           
  Partials        4        4           

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

[petrelharp]

Looking at the "details" on the failing tests (above), there's one failing test:

=========================== short test summary info ============================
FAILED tests/test_dfes.py::TestDroMelGamma_H17::test_mutation_types_match - AssertionError: assert False
 +  where False = <function allclose at 0x7ff425f22c30>([-0.000396, 0.33], [-0.000792, 0.33])
 +    where <function allclose at 0x7ff425f22c30> = np.allclose
 +    and   [-0.000396, 0.33] = MutationType(dominance_coeff=0.5, distribution_type='g', distribution_args=[-0.000396, 0.33], convert_to_substitution=True, dominance_coeff_list=None, dominance_coeff_breaks=None).distribution_args
 +    and   [-0.000792, 0.33] = MutationType(dominance_coeff=0.5, distribution_type='g', distribution_args=[-0.000792, 0.33], convert_to_substitution=True, dominance_coeff_list=None, dominance_coeff_breaks=None).distribution_args
===== 1 failed, 2599 passed, 53 skipped, 20 warnings in 349.21s (0:05:49) ======

You can run (only) this test locally by doing

python -m pytest tests/test_dfes.py::TestDroMelGamma_H17::test_mutation_types_match

• it should give the same error on your computer. And, the discrepancy looks like a factor of two! Could you double-check, and if you still agree with your version, briefly explain?

[clararehmann]

I get an error locally, I think it is a factor of two thing:

In my code, I use a -2 multiplier for the gamma_mean parameter to adjust it for SLiM, since per my reading it seems like the authors are using DaDi for their simulations (Text S2):

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],
    ...
    )

And in the original code for the DFE the mean is unscaled (though otherwise numerically the same):

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],
    )

I'm not sure which is more correct based on the paper!

[petrelharp]

I'm not sure which is more correct based on the paper!

Hm, why aren't you sure? If they're using dadi, there should be a factor of 2 (as you say); is it not clear that the parameters are from dadi? (I'm not looking at Text S2 right now; you could paste in a screenshot of the paragraph in question?)

[clararehmann]

This is the main text (referring to Text S2) where they say they use DaDi to infer s:

Based on this, it seems to me that there should be a factor of two in our code to convert to SLiM units?

[petrelharp]

Ah, I see - they used dadi to fit the demographic model, but then they "infer the DFE from the nonsynonymous SFS". Looks like we need to have a look in that Appendix (and, hopefully not the code).

[petrelharp]

Ah-ha: but in Text S2, they say

... in other words, they used dadi to do the calculations having to do with selection also. So, I agree with you about the factor of two.

[petrelharp]

So - could you fix up the original implementation with this factor of two as well? (And, the way you've coded it is better: writing down explicitly the shape/scale values that are in the paper and multiplying them to get the mean.) Then this should pass!

[clararehmann]

lol - I just took the same screenshot! I'll fix the implementation in the catalog as well :)

[petrelharp]

hah I was faster =)

[petrelharp]

Yay! Have a look at my suggestion; looks like we can merge once you accept that (or not?).