Merge pull request #1051 from andrewkern/contig_neut_test

added a function to test if a contig is neutral

stdpopsim/engines.py

@@ -237,6 +237,12 @@ def simulate(
                 del kwargs["random_seed"]
             else:
                 raise ValueError("Cannot set both seed and random_seed")
+        # test to make sure contig is fully neutral
+        if not contig.is_neutral():
+            raise ValueError(
+                "Contig had non neutral mutation types "
+                "but you are using the msprime engine"
+            )
 
         # TODO: remove this after a release or two. See #745.
         self._warn_zigzag(demographic_model)

stdpopsim/ext/selection.py

@@ -95,6 +95,12 @@ def __attrs_post_init__(self):
             "s": [],  # script types should just printout arguments
         }[self.distribution_type]
 
+    def is_neutral(self):
+        """
+        tests if the mutation type is strictly neutral
+        """
+        return self.distribution_type == "f" and self.distribution_args[0] == 0
+
 
 @attr.s
 class GenerationAfter(object):

stdpopsim/genomes.py

@@ -436,6 +436,13 @@ def fully_neutral(self, slim_mutations=False, convert_to_substitution=True):
             proportions=[1 if slim_mutations else 0],
         )
 
+    def is_neutral(self):
+        """
+        returns true if the contig has no non-neutral mutation
+        types
+        """
+        return all(mt.is_neutral() for mt in self.mutation_types)
+
     def __str__(self):
         gmap = "None" if self.genetic_map is None else self.genetic_map.id
         s = (

tests/test_engines.py

@@ -4,6 +4,7 @@
 import stdpopsim
 import msprime
 import pytest
+import numpy as np
 
 
 class TestEngineAPI:
@@ -104,3 +105,41 @@ def test_msprime_seed(self):
         sim_seed = engine.simulate(model, contig, samples, seed=1)
         sim_random_seed = engine.simulate(model, contig, samples, random_seed=1)
         assert sim_seed.tables.edges == sim_random_seed.tables.edges
+
+    def test_non_neutral_contig(self):
+        species = stdpopsim.get_species("HomSap")
+        model = species.get_demographic_model("AshkSub_7G19")
+        samples = model.get_samples(10)
+        contig = stdpopsim.Contig.basic_contig(length=100)
+        contig.clear_genomic_mutation_types()
+        props = [1]
+        mt = [stdpopsim.ext.MutationType(distribution_type="f", distribution_args=[1])]
+        dfes = [
+            stdpopsim.DFE(
+                id=str(0),
+                description="test",
+                long_description="test test",
+                proportions=props,
+                mutation_types=mt,
+            )
+        ]
+        contig.add_DFE(intervals=np.array([[0, 50]]), DFE=dfes[0])
+        engine = stdpopsim.get_engine("msprime")
+        with pytest.raises(ValueError):
+            engine.simulate(model, contig, samples, seed=1)
+        # okay now change selection coefficient to neutral
+        contig.clear_genomic_mutation_types()
+        props = [1]
+        mt = [stdpopsim.ext.MutationType(distribution_type="f", distribution_args=[0])]
+        dfes = [
+            stdpopsim.DFE(
+                id=str(0),
+                description="test",
+                long_description="test test",
+                proportions=props,
+                mutation_types=mt,
+            )
+        ]
+        contig.add_DFE(intervals=np.array([[0, 50]]), DFE=dfes[0])
+        engine = stdpopsim.get_engine("msprime")
+        engine.simulate(model, contig, samples, seed=1)

tests/test_genomes.py

@@ -111,6 +111,69 @@ def test_add_DFE(self):
         assert len(contig.genomic_element_types) == 2
         assert len(contig.mutation_types) == len(mt)
 
+    def test_is_neutral(self):
+        contig = stdpopsim.Contig.basic_contig(length=100)
+        contig.clear_genomic_mutation_types()
+        props = [0.3, 0.7]
+        mt = [
+            stdpopsim.ext.MutationType(distribution_type="f", distribution_args=[1])
+            for _ in props
+        ]
+        dfes = [
+            stdpopsim.DFE(
+                id=str(j),
+                description="test",
+                long_description="test test",
+                proportions=props,
+                mutation_types=mt,
+            )
+            for j in range(2)
+        ]
+        contig.add_DFE(intervals=np.array([[10, 30], [50, 100]]), DFE=dfes[0])
+        contig.add_DFE(intervals=np.array([[30, 40]]), DFE=dfes[1])
+        assert not contig.is_neutral()
+
+    def test_is_neutral2(self):
+        contig = stdpopsim.Contig.basic_contig(length=100)
+        contig.clear_genomic_mutation_types()
+        props = [0.3, 0.7]
+        mt = [stdpopsim.ext.MutationType() for _ in props]
+        dfes = [
+            stdpopsim.DFE(
+                id=str(j),
+                description="test",
+                long_description="test test",
+                proportions=props,
+                mutation_types=mt,
+            )
+            for j in range(2)
+        ]
+        contig.add_DFE(intervals=np.array([[10, 30], [50, 100]]), DFE=dfes[0])
+        contig.add_DFE(intervals=np.array([[30, 40]]), DFE=dfes[1])
+        assert contig.is_neutral()
+
+    def test_is_neutral3(self):
+        contig = stdpopsim.Contig.basic_contig(length=100)
+        contig.clear_genomic_mutation_types()
+        props = [0.3, 0.7]
+        mt = [
+            stdpopsim.ext.MutationType(distribution_type="e", distribution_args=[1])
+            for _ in props
+        ]
+        dfes = [
+            stdpopsim.DFE(
+                id=str(j),
+                description="test",
+                long_description="test test",
+                proportions=props,
+                mutation_types=mt,
+            )
+            for j in range(2)
+        ]
+        contig.add_DFE(intervals=np.array([[10, 30], [50, 100]]), DFE=dfes[0])
+        contig.add_DFE(intervals=np.array([[30, 40]]), DFE=dfes[1])
+        assert not contig.is_neutral()
+
 
 class TestAll_DFE_Models:
     """

tests/test_slim_engine.py

@@ -965,6 +965,16 @@ def test_bad_distribution_args_l(self):
                     distribution_type="l", distribution_args=distribution_args
                 )
 
+    def test_netural_mutation_type(self):
+        assert stdpopsim.ext.MutationType().is_neutral()
+
+    def test_not_neutral_mutation_type(self):
+        mt = stdpopsim.ext.MutationType(
+            distribution_type="f",
+            distribution_args=[1],
+        )
+        assert not mt.is_neutral()
+
 
 @pytest.mark.skipif(IS_WINDOWS, reason="SLiM not available on windows")
 class TestDrawMutation(PiecewiseConstantSizeMixin):