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

Issue:
change BosTau/HolsteinFriesian_1M13 to be piecewise constant
See #579 for background.
write parameter table for BosTau/HolsteinFriesian_1M13
The BosTau/HolsteinFriesian_1M13 parameter table for the catalog is missing.
https://stdpopsim.readthedocs.io/en/latest/development.html#write-parameter-table

Needs #601 to be resolved first.


Discussion / hints:

@grahamgower Do you want to assign this to me?

Fix patch:
diff --git a/docs/parameter_tables/BosTau/HolsteinFriesian_1M13.csv b/docs/parameter_tables/BosTau/HolsteinFriesian_1M13.csv
new file mode 100644
index 000000000..fa1ca4c9a
--- /dev/null
+++ b/docs/parameter_tables/BosTau/HolsteinFriesian_1M13.csv
@@ -0,0 +1,26 @@
+Population size,"62,000",Ancestral population size
+Population size,"17,000",Pop. size during 1st time interval
+Population size,"10,000",Pop. size during 2nd time interval
+Population size,"7,000",Pop. size during 3rd time interval
+Population size,"3,500",Pop. size during 4th time interval
+Population size,"2,500",Pop. size during 5th time interval
+Population size,"2,000",Pop. size during 6th time interval
+Population size,"1,500",Pop. size during 7th time interval
+Population size,"1,000",Pop. size during 8th time interval
+Population size,"350",Pop. size during 9th time interval
+Population size,"250",Pop. size during 10th time interval
+Population size,"120",Pop. size during 11th time interval
+Population size,"90",Pop. size during 12th time interval
+Time (gen.),"33,154",Begining of 1st time interval
+Time (gen.),"3,354",Begining of 2nd time interval
+Time (gen.),"2,354",Begining of 3rd time interval
+Time (gen.),"1,754",Begining of 4th time interval
+Time (gen.),"654",Begining of 5th time interval
+Time (gen.),"454",Begining of 6th time interval
+Time (gen.),"154",Begining of 7th time interval
+Time (gen.),"24",Begining of 8th time interval
+Time (gen.),"18",Begining of 9th time interval
+Time (gen.),"12",Begining of 10th time interval
+Time (gen.),"6",Begining of 11th time interval
+Time (gen.),"3",Begining of 12th time interval
+Generation time (yrs.),5,Generation time
diff --git a/stdpopsim/catalog/BosTau/__init__.py b/stdpopsim/catalog/BosTau/__init__.py
index 48461dd7f..51fc1bbbd 100644
--- a/stdpopsim/catalog/BosTau/__init__.py
+++ b/stdpopsim/catalog/BosTau/__init__.py
@@ -1,5 +1,5 @@
 """
-Catalog definitions for Bos Taurus
+Catalog definitions for bos_taurus
 """
 from . import species  # noqa: F401
 from . import demographic_models  # noqa: F401
diff --git a/stdpopsim/catalog/BosTau/demographic_models.py b/stdpopsim/catalog/BosTau/demographic_models.py
index e42b761ed..9129a027a 100644
--- a/stdpopsim/catalog/BosTau/demographic_models.py
+++ b/stdpopsim/catalog/BosTau/demographic_models.py
@@ -1,31 +1,28 @@
 import msprime
-
 import stdpopsim
 
 _species = stdpopsim.get_species("BosTau")
 
-# Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence,
-# with Correction for Sequence Errors.
-_MacLeodEtAl = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/molbev/mst125",
-    year="2013",
-    author="MacLeod et al.",
-)
-
 
-def _HolsteinFriesan_1M13():
+def _HolsteinFriesian_1M13():
     id = "HolsteinFriesian_1M13"
     description = "Piecewise size changes in Holstein-Friesian cattle."
     long_description = """
     The piecewise-constant population size model of Holstein-Friesian cattle
-    from MacLeod et al. 2013. Population sizes were estimated from inferred
-    runs of homozygosity, with parameter values taken from Figure 4A by visual
-    inspection of the plots.
+    from MacLeod et al. (2013). Population sizes were estimated from inferred
+    runs of homozygosity, with parameter values taken from Figure 4A and Table S1.
     """
     populations = [
-        stdpopsim.Population(id="Holstein_Friesian", description="Holstein Friesian"),
+        stdpopsim.Population(id="Holstein_Friesian", description="Holstein-Friesian"),
+    ]
+    citations = [
+        stdpopsim.Citation(
+            author="MacLeod et al.",
+            year=2013,
+            doi="https://doi.org/10.1093/molbev/mst125",
+            reasons={stdpopsim.CiteReason.DEM_MODEL},
+        )
     ]
-    citations = [_MacLeodEtAl.because(stdpopsim.CiteReason.DEM_MODEL)]
 
     return stdpopsim.DemographicModel(
         id=id,
@@ -35,64 +32,64 @@ def _HolsteinFriesan_1M13():
         citations=citations,
         generation_time=_species.generation_time,
         population_configurations=[
+            #      3 generations at    90,     1-     3
             msprime.PopulationConfiguration(
-                initial_size=90, growth_rate=0.0166, metadata=populations[0].asdict()
+                initial_size=90, metadata=populations[0].asdict()
             )
         ],
         # Here 'time' should be in generation notation ie. how many
         # generations ago were that Ne (effective population size)
-        # and growth rate.
-        # Growth rate is "per generation exponential growth rate":
-        # -alpha= [ln(initial_pop_size/next_stage_pop_size)/generation_span_in_years]
-        # For example: ln(90/120)/3= -0.095894024
         demographic_events=[
+            #      3 generations at   120,     4-     6
             msprime.PopulationParametersChange(
-                time=1,
-                initial_size=90,
-                growth_rate=-0.095894024,
-                population_id=0,
-            ),  # Ne 90 to 120
-            msprime.PopulationParametersChange(
-                time=4, growth_rate=-0.24465639, population_id=0
-            ),  # Ne 120 to 250
-            msprime.PopulationParametersChange(
-                time=7, growth_rate=-0.0560787, population_id=0
-            ),  # Ne 250 to 350
+                time=4, initial_size=120, population_id=0
+            ),
+            #      6 generations at   250,     7-    12
             msprime.PopulationParametersChange(
-                time=13, growth_rate=-0.1749704, population_id=0
-            ),  # Ne 350 to 1000
+                time=7, initial_size=250, population_id=0
+            ),
+            #      6 generations at   350,    13-    18
             msprime.PopulationParametersChange(
-                time=19, growth_rate=-0.0675775, population_id=0
-            ),  # Ne 1000 to 1500
+                time=13, initial_size=350, population_id=0
+            ),
+            #      6 generations at  1000,    19-    24
             msprime.PopulationParametersChange(
-                time=25, growth_rate=-0.0022129, population_id=0
-            ),  # Ne 1500 to 2000
+                time=19, initial_size=1000, population_id=0
+            ),
+            #    130 generations at  1500,    25-   154
             msprime.PopulationParametersChange(
-                time=155, growth_rate=-0.0007438, population_id=0
-            ),  # Ne 2000 to 2500
+                time=25, initial_size=1500, population_id=0
+            ),
+            #    200 generations at  2000,   155-   454
             msprime.PopulationParametersChange(
-                time=455, growth_rate=-0.0016824, population_id=0
-            ),  # Ne 2500 to 3500
+                time=155, initial_size=2000, population_id=0
+            ),
+            #    200 generations at  2500,   455-   654
             msprime.PopulationParametersChange(
-                time=655, growth_rate=-0.0006301, population_id=0
-            ),  # Ne 3500 to 7000
+                time=455, initial_size=2500, population_id=0
+            ),
+            #   1100 generations at  3500,   655-  1754
             msprime.PopulationParametersChange(
-                time=1755, growth_rate=-0.0005945, population_id=0
-            ),  # Ne 7000 to 10000
+                time=655, initial_size=3500, population_id=0
+            ),
+            #    600 generations at  7000,  1755-  2354
             msprime.PopulationParametersChange(
-                time=2355, growth_rate=-0.0005306, population_id=0
-            ),  # Ne 10000 to 17000
+                time=1755, initial_size=7000, population_id=0
+            ),
+            #   1000 generations at 10000,  2355-  3354
             msprime.PopulationParametersChange(
-                time=3355, growth_rate=-0.0000434, population_id=0
-            ),  # Ne 17000 to 62000
+                time=2355, initial_size=10000, population_id=0
+            ),
+            #  29800 generations at 17000,  3355- 33154
             msprime.PopulationParametersChange(
-                time=33155, growth_rate=-0.0000, population_id=0
-            ),  # Ne 62000 (model has "coalesced")
+                time=3355, initial_size=17000, population_id=0
+            ),
+            # 900000 generations at 62000, 33155-933154
             msprime.PopulationParametersChange(
-                time=933155, growth_rate=-0.0, population_id=0
+                time=33155, initial_size=62000, population_id=0
             ),
         ],
     )
 
 
-_species.add_demographic_model(_HolsteinFriesan_1M13())
+_species.add_demographic_model(_HolsteinFriesian_1M13())
diff --git a/stdpopsim/catalog/BosTau/species.py b/stdpopsim/catalog/BosTau/species.py
index 85270c2c9..386f46b2e 100644
--- a/stdpopsim/catalog/BosTau/species.py
+++ b/stdpopsim/catalog/BosTau/species.py
@@ -1,39 +1,44 @@
 import collections
-
 import stdpopsim
+
 from . import genome_data
 
 # De novo assembly of the cattle reference genome with single-molecule sequencing.
 _RosenEtAl = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/gigascience/giaa021",
-    year="2020",
     author="Rosen et al.",
+    year=2020,
+    doi="https://doi.org/10.1093/gigascience/giaa021",
+    reasons={stdpopsim.CiteReason.ASSEMBLY},
 )
 
 # Frequency of mosaicism points towards mutation-prone early cleavage
 # cell divisions in cattle.
 _HarlandEtAl = stdpopsim.Citation(
     author="Harland et al.",
-    year="2017",
-    # BioRxiv preprint
+    year=2017,
     doi="https://doi.org/10.1101/079863",
+    reasons={stdpopsim.CiteReason.MUT_RATE},
 )
 
 # Cattle Sex-Specific Recombination and Genetic Control from a
 # Large Pedigree Analysis.
 _MaEtAl = stdpopsim.Citation(
     author="Ma et al.",
-    year="2015",
+    year=2015,
     doi="https://doi.org/10.1371/journal.pgen.1005387",
+    reasons={stdpopsim.CiteReason.REC_RATE},
 )
 
 # Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence,
 # with Correction for Sequence Errors.
 _MacLeodEtAl = stdpopsim.Citation(
-    doi="https://doi.org/10.1093/molbev/mst125",
-    year="2013",
     author="MacLeod et al.",
-    reasons={stdpopsim.CiteReason.GEN_TIME, stdpopsim.CiteReason.POP_SIZE},
+    year=2013,
+    doi="https://doi.org/10.1093/molbev/mst125",
+    reasons={
+        stdpopsim.CiteReason.GEN_TIME,
+        stdpopsim.CiteReason.POP_SIZE,
+    },
 )
 
 # Recombination rate has been derived from dairy cattle crossovers
@@ -68,20 +73,20 @@
 _genome = stdpopsim.Genome(
     chromosomes=_chromosomes,
     citations=[
-        _HarlandEtAl.because(stdpopsim.CiteReason.MUT_RATE),
-        _MaEtAl.because(stdpopsim.CiteReason.REC_RATE),
-        _RosenEtAl.because(stdpopsim.CiteReason.ASSEMBLY),
+        _RosenEtAl,
+        _HarlandEtAl,
+        _MaEtAl,
     ],
 )
 
 _species = stdpopsim.Species(
     id="BosTau",
     ensembl_id="bos_taurus",
-    name="Bos Taurus",
+    name="Bos taurus",
     common_name="Cattle",
     genome=_genome,
     generation_time=5,
-    population_size=62000,
+    population_size=90,  # most recent Ne in _MacLeodEtAl
     citations=[_MacLeodEtAl],
 )
 


Test patch:
diff --git a/tests/test_BosTau.py b/tests/test_BosTau.py
new file mode 100644
index 000000000..fe6e58d9f
--- /dev/null
+++ b/tests/test_BosTau.py
@@ -0,0 +1,113 @@
+import pytest
+
+import stdpopsim
+from tests import test_species
+
+
+class TestSpecies(test_species.SpeciesTestBase):
+
+    species = stdpopsim.get_species("BosTau")
+
+    def test_ensembl_id(self):
+        assert self.species.ensembl_id == "bos_taurus"
+
+    def test_name(self):
+        assert self.species.name == "Bos taurus"
+
+    def test_common_name(self):
+        assert self.species.common_name == "Cattle"
+
+    def test_qc_population_size(self):
+        assert self.species.population_size == 90  # most recent Ne in _MacLeodEtAl
+
+    def test_qc_generation_time(self):
+        assert self.species.generation_time == 5
+
+
+class TestGenome(test_species.GenomeTestBase):
+
+    genome = stdpopsim.get_species("BosTau").genome
+
+    def test_basic_attributes(self):
+        nchrom = 31  # 29 + X + MT
+        assert len(self.genome.chromosomes) == nchrom
+
+    @pytest.mark.skip("Recombination rate QC not done yet")
+    @pytest.mark.parametrize(
+        ["name", "rate"],
+        {
+            "1": -1,
+            "2": -1,
+            "3": -1,
+            "4": -1,
+            "5": -1,
+            "6": -1,
+            "7": -1,
+            "8": -1,
+            "9": -1,
+            "10": -1,
+            "11": -1,
+            "12": -1,
+            "13": -1,
+            "14": -1,
+            "15": -1,
+            "16": -1,
+            "17": -1,
+            "18": -1,
+            "19": -1,
+            "20": -1,
+            "21": -1,
+            "22": -1,
+            "23": -1,
+            "24": -1,
+            "25": -1,
+            "26": -1,
+            "27": -1,
+            "28": -1,
+            "29": -1,
+            "X": -1,
+            "MT": -1,
+        }.items(),
+    )
+    def test_recombination_rate(self, name, rate):
+        assert pytest.approx(rate, self.genome.get_chromosome(name).recombination_rate)
+
+    @pytest.mark.skip("Mutation rate QC not done yet")
+    @pytest.mark.parametrize(
+        ["name", "rate"],
+        {
+            "1": -1,
+            "2": -1,
+            "3": -1,
+            "4": -1,
+            "5": -1,
+            "6": -1,
+            "7": -1,
+            "8": -1,
+            "9": -1,
+            "10": -1,
+            "11": -1,
+            "12": -1,
+            "13": -1,
+            "14": -1,
+            "15": -1,
+            "16": -1,
+            "17": -1,
+            "18": -1,
+            "19": -1,
+            "20": -1,
+            "21": -1,
+            "22": -1,
+            "23": -1,
+            "24": -1,
+            "25": -1,
+            "26": -1,
+            "27": -1,
+            "28": -1,
+            "29": -1,
+            "X": -1,
+            "MT": 0,
+        }.items(),
+    )
+    def test_mutation_rate(self, name, rate):
+        assert pytest.approx(rate, self.genome.get_chromosome(name).mutation_rate)