Msprime: Add MsmsSweeps test to verification.py

[jeromekelleher]

As an additional check on our sweeps code, we'd like to verify against msms. Initially we'll use this to check our genic sweeps code. This will follow the same general approach as the current DiscoalSweeps class.

To make this work we need to add msms to the Makefile in the data directory. This should download and msms.

Once this basic infrastructure works, we should at ~10 examples of msms command lines which try out various parameter regimes that we support.

cc @andrewkern

[jeromekelleher]

cc @gbinux

[gbinux]

There are two modes of selection in msms: -SI goes backward allowing us to set end
frequency at specified ending time; -SF goes forward allowing us to set start frequency at
specified starting time.

Take -SI for example. The relevant options are -SI start_time 1 start_freq -SAA sAA -SaA saA -Sp sel_pos. When converting them to msprime options, I found some of them
might not have the counterparts:

1. msms could specify selection strengths for both homozygous (sAA) for
   selected allele and heterozygous (saA). In msprime/discoal, how could we calculate/specify the saA or h parameters?
2. msms allows us to set the time when selective sweep starts or ends ( not
   both). msprime/ancestry.py: class SweepGenicSelection likely does not have a
   parameter to specify the selection start/end time? Is this correct? @andrewkern

[andrewkern]

Take -SI for example. The relevant options are -SI start_time 1 start_freq -SAA sAA -SaA saA -Sp sel_pos. When converting them to msprime options, I found some of them
might not have the counterparts:

1. msms could specify selection strengths for both homozygous (sAA) for
   selected allele and heterozygous (saA). In msprime/discoal, how could we calculate/specify the saA or h parameters?

right now we just have genic selection implemented so h=0.5

1. msms allows us to set the time when selective sweep starts or ends ( not
   both). msprime/ancestry.py: class SweepGenicSelection likely does not have a
   parameter to specify the selection start/end time? Is this correct? @andrewkern

yes we have implemented the time when it ends (looking forwards in time). there are a few examples of that in test_models.py

[gbinux]

Got it. Time parameter is not within the model class. Thanks a lot! @andrewkern

[gbinux]

There seems to be some inconsitency between the hudson model and the SweepGenicSelection model in terms of the scale factor of mutation and recombination rates. When hudson model is specified in msprime.simulate method, mutation rate and recombination need to be provided as raw numbers (per generation per base); while with SweepGenicSelection model, mutation rate and recombination rate seem to be scaled by 4 * N (hint from verification.py class DiscoalTest method _discoal_str_to_simulation). In one test case, I found SweepGenicSelection model with rates scaled by 4N generates more segsites hudson with raw rates; if usingSweepGenicSelection` model with rates scaled by 2N, then two model generated comparable numer of segsites. I will do some more exploration about this later. @andrewkern

[andrewkern]

so yes there is an issue with scaling between the models, but it's with respect to time rather than mutation and recombination rates. if you look in the msprime model, there is currently a factor of 4N being used

https://github.com/tskit-dev/msprime/blob/dab983c4bebb188d10ce6ff2b2f07b0636d25fc0/lib/msprime.c#L3639-L3645

the sweep implementation I have ported from discoal originally used a 2N time scale internally, and while I thought I had taken care of that it didn't work. i _think_ we should be able to multiply by a scalar here:

https://github.com/tskit-dev/msprime/blob/dab983c4bebb188d10ce6ff2b2f07b0636d25fc0/lib/msprime.c#L4612

however i'm struggling to get msprime and discoal to match

[andrewkern]

Looks good with msms 5 300 -t 200 -r 200 500000 -N 10000 -SaA 5000 -SAA 10000 -N 10000 simulation if use 2N scale for SweepGenicSelection model.

huh... what do you mean by "use 2N scale"?

[gbinux]

oops, made a mistake...

[gbinux]

Ne=0.25 ( instead of Ne=refsize) seems critical for msprime to generate stats comparable to those from msms but I am not sure why we need to useNe=0.25 but not Ne=refsize @andrewkern

https://github.com/tskit-dev/msprime/blob/66972a7d213a0625d5dc3cdb0f370726dc924755/verification.py#L746

[andrewkern]

Yes setting Ne=0.25 allows the results to match canonical coalescent scaling. You'll see the same thing with the ms verification (via mspms cli).

[molpopgen]

Does that mean the mutation/recombination rates are 4Nu and 4Nr for the sweep models? If so, would it be better to change them to match the standard msprime models?

[andrewkern]

Well not really cause the time scaling for sweeps is still screwed.... @molpopgen got a chance to chat with me about this tomorrow? I could use talking it over with someone.

Basic issue is that sweeps are happening in 2N time scaling, but want it to match the 4N currently being used by msprime and I'm running into weirdness

[molpopgen]

Yeah I could probably find some time to chat tomorrow. Send me a message on slack and we'll figure it out.

[molpopgen]

My days are tending to start a little later in the morning, just FYI

[gbinux]

1

My new MsmsSweeps test class started to work "normally" with when setting recombination_rate=rho/(num_sites-1), mutation_rate=theta/num_sites and Ne=0.25 for the sweep model.

With relative smaller small sizes, the stats between msms and msprime simulations were similar in both hudson and geneic selection models. I like this! @andrewkern

2

But when I increased sample size, say 300, the stats for msms and msprime simulations became very different in the genic selection model ( still similar in Hudson model).

I am thinking that internally there are some big differences between msms and discoal/msprime in the methods of generating frequency trajectory. So it is reasonable to see very different stats between them in some cases? @jeromekelleher @andrewkern

3

There seems to be a bug. I saw the following error with some particular parameters at some time but the same parameters could run normally at a different time. And an interesting observation is that the error only happened when simulating selective sweep?! I am not sure whether I should put this into a new issue.

    python3: lib/msprime.c:5845: msp_std_common_ancestor_event: Assertion `y_node != NULL' failed.
    Aborted

4

Currently, I just tested sweeps in 1 population with constant Ne. I saw multiple sweeps in test_models.py. I can try this next.
@andrewkern May I have your suggestions for other features to be included in MsmsSweeps test class? And how can I push my code based on the version of verification.py and data/Makefile in pull request #1104?

[andrewkern]

I think the best thing to do is wait for #1104 to be merged once i get the scaling, and then you should start a PR. this code is still a work in progress.

btw can you post the sweep comparison figure a bit larger? curious to see what you are getting out

[gbinux]

btw can you post the sweep comparison figure a bit larger? curious to see what you are getting out

I just uploaded larger pictures and updated the original post above. @andrewkern

I think the best thing to do is wait for #1104 to be merged once i get the scaling, and then you should start a PR. this code is still a work in progress.

Agreed. I will wait for #1104 to be merged then. Look forward to your updates. :+1:

[jeromekelleher]

This looks like great progress @gbinux! I think we should merge #1104 ASAP (without trying to fix the underlying problems) and get a PR opened with the msms tests. The more tests we have the more information we'll have when trying to pin down what the underlying issues are.

[jeromekelleher]

Closed in #1109.