Repository: popsim-consortium/stdpopsim (language: python)

Issue:
HomSap/Zigzag_1S14 population sizes don't match Schiffels & Durbin (2014)
**Short description of the problem:**
The population sizes in the `Zigzag_1S14` demographic model swing between 1431 and 14312. However, the values in [Schiffels & Durbin (2014)](https://doi.org/10.1038/ng.3015) are five times higher.

**How to reproduce the problem:**
Compare the zigzag parameters in [stdpopsim's catalog](https://stdpopsim.readthedocs.io/en/latest/catalog.html#sec_catalog_homsap_models_zigzag_1s14) with Figure 2 in Schiffels & Durbin (2014). Also see the `ms` command in the supplementary material, section 7, pg. 25.

This appears to have been deliberate, according to #105, but is not documented, and doesn't seem to have come up during the model QC process (#498).

**Suggested fixes:**
The original `ms` command given in the supplement includes a `-eN 0 5` command, which has not been reproduced in the stdpopsim model. From the ms manual:
`-eN t x` `Set all subpop’s to size x ∗ N0 and growth rates to zero.`

We should at the very least document the difference between the stdpopsim ZigZag model, and that in the original publication.

CC: @Chris1221, @apragsdale.

[EDIT: [translation of the ms command, as stolen from the QC implementation](https://gist.github.com/grahamgower/7a2c0115493547e5ddc02887030abee9)]


Discussion / hints:
Yep, that looks like a mistake on our part. We should probably fix those values in our implementation to follow the original publication, instead of just documenting it. Maybe the next major release could emit a warning the first time someone uses the zigzag model, that the previous version of stdpopsim had sizes that were 5x smaller than intended?

I'm sorry for missing that in the QC. I'll be very happy the day that I never have to translate an ms command again.
Happy day indeed! It'll come: https://github.com/popsim-consortium/demes-python/pull/102

(more quickly, if someone would like to pick this PR up...)
I agree, that does look like a mistake in our implementation. My apologies for missing that initially, my ms-fu is pretty weak and I completely didn't see the scaling factor. Comparing it to the literature values would have been a good sanity check that I also missed... I agree with @apragsdale's warning suggestion, it's important that people know that it was incorrect!  
Also, (I think) it looks like this change may need to be propagated to the values in the Zigzag model description in the demes-docs gallery: https://popsim-consortium.github.io/demes-docs/main/gallery.html#zigzag-model
> Also, (I think) it looks like this change may need to be propagated to the values in the Zigzag model description in the demes-docs gallery: https://popsim-consortium.github.io/demes-docs/main/gallery.html#zigzag-model

Yep, thanks @Chris1221, I've opened an issues in Demes for this too: https://github.com/popsim-consortium/demes-python/issues/179.

> Also, (I think) it looks like this change may need to be propagated to the values in the Zigzag model description in the demes-docs gallery: https://popsim-consortium.github.io/demes-docs/main/gallery.html#zigzag-model

Yes, for sure. I think I used our implementation from stdpopsim to translate over to demes. Great example of how easily these errors can propagate.. similar to the OOA model issue from earlier this year. We'll need fix it there as well.
Great to catch this issue, but I don't think it's that serious in the scheme of things since it's a toy model. It's just another example of how easy it is to make mistakes with today's toolchain. So, no worries @Chris1221 and @apragsdale!

In terms of compatibility, I'm not sure what to do. I guess emitting a warning when it's used is probably the best we can do, which we remove after (say) 6 months.

Fix patch:
diff --git a/stdpopsim/catalog/HomSap/__init__.py b/stdpopsim/catalog/HomSap/__init__.py
index eccc390c5..042e45e51 100644
--- a/stdpopsim/catalog/HomSap/__init__.py
+++ b/stdpopsim/catalog/HomSap/__init__.py
@@ -895,7 +895,7 @@ def _zigzag():
     ]
 
     generation_time = 30
-    N0 = 14312
+    N0 = 5 * 14312
 
     g_1 = 0.02302578256
     t_1 = 33.333334976
diff --git a/stdpopsim/engines.py b/stdpopsim/engines.py
index a423a905a..bf461a296 100644
--- a/stdpopsim/engines.py
+++ b/stdpopsim/engines.py
@@ -1,4 +1,5 @@
 import logging
+import warnings
 
 import attr
 import msprime
@@ -101,6 +102,15 @@ def get_version(self):
         """
         raise NotImplementedError()
 
+    def _warn_zigzag(self, demographic_model):
+        if demographic_model.id == "Zigzag_1S14":
+            warnings.warn(
+                "In stdpopsim <= 0.1.2, the Zigzag_1S14 model produced population "
+                "sizes 5x lower than those in Schiffels & Durbin (2014). "
+                "The population sizes have now been corrected. For details see: "
+                "https://github.com/popsim-consortium/stdpopsim/issues/745"
+            )
+
 
 class _MsprimeEngine(Engine):
     id = "msprime"  #:
@@ -181,6 +191,9 @@ def simulate(
             else:
                 raise ValueError("Cannot set both seed and random_seed")
 
+        # TODO: remove this after a release or two. See #745.
+        self._warn_zigzag(demographic_model)
+
         ts = msprime.simulate(
             samples=samples,
             recombination_map=contig.recombination_map,
diff --git a/stdpopsim/qc/HomSap.py b/stdpopsim/qc/HomSap.py
index aedb0d706..6585a7904 100644
--- a/stdpopsim/qc/HomSap.py
+++ b/stdpopsim/qc/HomSap.py
@@ -1103,11 +1103,12 @@ def ZigZag():
     mu = 1.25e-8
     L = 10000000
     theta = 7156
-    Ne = theta / 4 / mu / L
+    Ne = theta / 4 / mu / L  # ancestral size for rescaling times/rates
+    N0 = 5 * Ne  # multiply pop sizes to match 'ms -eN 0 5'
 
     population_configurations = [
         msprime.PopulationConfiguration(
-            initial_size=Ne, growth_rate=0, metadata={"sampling_time": 0}
+            initial_size=N0, growth_rate=0, metadata={"sampling_time": 0}
         )
     ]
 
@@ -1133,7 +1134,7 @@ def ZigZag():
             time=0.596236 * 4 * Ne,
             growth_rate=0,
             population_id=0,
-            initial_size=0.1 * Ne,
+            initial_size=0.1 * N0,
         )
     )
 
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index 6a268c864..b7bde9800 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -1070,6 +1070,9 @@ def simulate(
                 )
             )
 
+        # TODO: remove this after a release or two. See #745.
+        self._warn_zigzag(demographic_model)
+
         run_slim = not slim_script
 
         # Ensure only "weighted" mutations are introduced by SLiM.


Test patch:
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 37dc1e9c8..e8d3a980d 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -244,10 +244,6 @@ def test_ragsdale_archaic(self):
         cmd = "HomSap -c chr1 -l0.01 -d OutOfAfricaArchaicAdmixture_5R19 10"
         self.verify(cmd, num_samples=10)
 
-    def test_schiffels_zigzag(self):
-        cmd = "HomSap -c chr1 -l0.01 -d Zigzag_1S14 2"
-        self.verify(cmd, num_samples=2)
-
     def test_dromel_constant(self):
         cmd = "DroMel -c 2L -l0.001 4"
         self.verify(cmd, num_samples=4)
diff --git a/tests/test_models.py b/tests/test_models.py
index 21ffa1f4a..032f8b180 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -606,3 +606,14 @@ def test_population_construction_no_popconfig(self):
             id="A", description="A", long_description="A", generation_time=1
         )
         self.assertEqual(dm.populations, [])
+
+
+class TestZigZagWarning(unittest.TestCase):
+    def test_zigzag_produces_warning(self):
+        species = stdpopsim.get_species("HomSap")
+        model = species.get_demographic_model("Zigzag_1S14")
+        contig = species.get_contig("chr22")
+        samples = model.get_samples(10)
+        for engine in stdpopsim.all_engines():
+            with self.assertWarnsRegex(UserWarning, "Zigzag_1S14"):
+                engine.simulate(model, contig, samples, dry_run=True)