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

Issue:
add vaquita
As discussed today, vaquita has a good assembly and various demographic models. It even has a DFE inferred, but @klohmueller says that it doesn't fit significantly better than the human-inferred one added in #1505.

Here's the info:  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9881057/
and see Table S2 for demographic parameter estimates.




Discussion / hints:
@igronau has volunteered to do this! <3 <3

Fix patch:
diff --git a/README.md b/README.md
index 8fbfd63ec..221b7edd4 100644
--- a/README.md
+++ b/README.md
@@ -11,6 +11,10 @@ for in-development docs, see [this link](https://popsim-consortium.github.io/std
 
 Click [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/popsim-consortium/stdpopsim/main?filepath=stdpopsim_example.ipynb) to start an interactive Jupyter Notebook and start playing with `stdpopsim` now!
 
-If you use ``stdpopsim`` in your work, please cite our papers ([original](https://elifesciences.org/articles/54967) and [followup](https://www.biorxiv.org/content/10.1101/2022.10.29.514266v1)).
+If you use ``stdpopsim`` in your work, please cite our papers:
+
+* [Adrion, et al. (2020)](https://elifesciences.org/articles/54967)
+* [Lauterbur, et al. (2023)](https://elifesciences.org/reviewed-preprints/84874).
+
 See [here](https://popsim-consortium.github.io/stdpopsim-docs/stable/introduction.html#citations)
 for full citation details.
diff --git a/docs/introduction.rst b/docs/introduction.rst
index 861a0c4f7..15a4a8d09 100644
--- a/docs/introduction.rst
+++ b/docs/introduction.rst
@@ -77,20 +77,20 @@ major additions to the catalog:
     Travis J Struck, Simon Gravel, Ryan N Gutenkunst, Kirk E Lohmeuller,
     Peter L Ralph, Daniel R Schrider, Adam Siepel, Jerome Kelleher, Andrew D Kern (2020),
     *A community-maintained standard library of population genetic models*,
-    eLife 2020;9:e54967; doi: https://doi.org/10.7554/eLife.54967
-
-  - M. Elise Lauterbur, Mariz Izabel A. Cavassim, Ariella L. Gladstein, Graham
-    Gower, Nathaniel S. Pope, Georgia Tsambos, Jeffrey Adrion, Saurabh Belsare,
-    Arjun Biddanda, Victoria Caudill, Jean Cury, Ignacio Echevarria, Benjamin C.
-    Haller, Ahmed R. Hasan, Xin Huang, Leonardo Nicola Martin Iasi, Ekaterina
-    Noskova, Jana Obšteter, Vitor Antonio Corrêa Pavinato, Alice Pearson, David
-    Peede, Manolo F. Perez, Murillo F. Rodrigues, Chris C. R. Smith, Jeffrey P.
-    Spence, Anastasia Teterina, Silas Tittes, Per Unneberg, Juan Manuel Vazquez,
-    Ryan K. Waples, Anthony Wilder Wohns, Yan Wong, Franz Baumdicker, Reed A.
-    Cartwright, Gregor Gorjanc, Ryan N. Gutenkunst, Jerome Kelleher, Andrew D.
-    Kern, Aaron P. Ragsdale, Peter L. Ralph, Daniel R. Schrider, Ilan Gronau (2022).
+    eLife 9:e54967; doi: https://doi.org/10.7554/eLife.54967
+
+  - M. Elise Lauterbur, Maria Izabel A. Cavassim, Ariella L. Gladstein, Graham Gower,
+    Nathaniel S. Pope, Georgia Tsambos, Jeff Adrion, Saurabh Belsare, Arjun Biddanda,
+    Victoria Caudill, Jean Cury, Ignacio Echevarria, Benjamin C. Haller, Ahmed R. Hasan,
+    Xin Huang, Leonardo Nicola Martin Iasi, Ekaterina Noskova, Jana Obšteter,
+    Vitor Antonio Corrêa Pavinato, Alice Pearson, David Peede, Manolo F. Perez,
+    Murillo F. Rodrigues, Chris C. R. Smith, Jeffrey P. Spence, Anastasia Teterina,
+    Silas Tittes, Per Unneberg, Juan Manuel Vazquez, Ryan K. Waples, Anthony Wilder Wohns,
+    Yan Wong, Franz Baumdicker, Reed A. Cartwright, Gregor Gorjanc, Ryan N. Gutenkunst,
+    Jerome Kelleher, Andrew D. Kern, Aaron P. Ragsdale, Peter L. Ralph, Daniel R. Schrider,
+    Ilan Gronau (2023)
     *Expanding the stdpopsim species catalog, and lessons learned for realistic genome simulations*,
-    bioRxiv 2022.10.29.514266; doi: https://doi.org/10.1101/2022.10.29.514266
+    eLife 12:RP84874; doi: https://doi.org/10.7554/eLife.84874.2
 
 Bibtex records::
 
@@ -121,30 +121,23 @@ Bibtex records::
         publisher = {eLife Sciences Publications, Ltd},
     }
 
-    @article {Lauterbur2022.10.29.514266,
-        author = {Lauterbur, M. Elise and Cavassim, Mariz Izabel A. and
-            Gladstein, Ariella L. and Gower, Graham and Pope, Nathaniel S. and
-            Tsambos, Georgia and Adrion, Jeffrey and Belsare, Saurabh and Biddanda,
-            Arjun and Caudill, Victoria and Cury, Jean and Echevarria, Ignacio and
-            Haller, Benjamin C. and Hasan, Ahmed R. and Huang, Xin and Iasi,
-            Leonardo Nicola Martin and Noskova, Ekaterina and Ob{\v s}teter, Jana
-            and Pavinato, Vitor Antonio Corr{\^e}a and Pearson, Alice and Peede,
-            David and Perez, Manolo F. and Rodrigues, Murillo F. and Smith, Chris
-            C. R. and Spence, Jeffrey P. and Teterina, Anastasia and Tittes, Silas
-            and Unneberg, Per and Vazquez, Juan Manuel and Waples, Ryan K. and
-            Wohns, Anthony Wilder and Wong, Yan and Baumdicker, Franz and
-            Cartwright, Reed A. and Gorjanc, Gregor and Gutenkunst, Ryan N. and
-            Kelleher, Jerome and Kern, Andrew D. and Ragsdale, Aaron P. and Ralph,
-            Peter L. and Schrider, Daniel R. and Gronau, Ilan},
-        title = {Expanding the stdpopsim species catalog, and lessons learned
-            for realistic genome simulations},
-        elocation-id = {2022.10.29.514266},
-        year = {2022},
-        doi = {10.1101/2022.10.29.514266},
-        publisher = {Cold Spring Harbor Laboratory},
-        URL = {https://www.biorxiv.org/content/early/2022/10/31/2022.10.29.514266},
-        eprint = {https://www.biorxiv.org/content/early/2022/10/31/2022.10.29.514266.full.pdf},
-        journal = {bioRxiv}
+    @article{10.7554/eLife.84874.2,
+        title={Expanding the stdpopsim species catalog, and lessons learned for realistic genome simulations},
+	author = {M. Elise Lauterbur and Maria Izabel A. Cavassim and Ariella L. Gladstein and Graham Gower and
+           Nathaniel S. Pope and Georgia Tsambos and Jeff Adrion and Saurabh Belsare and Arjun Biddanda and
+           Victoria Caudill and Jean Cury and Ignacio Echevarria and Benjamin C. Haller and Ahmed R. Hasan and
+           Xin Huang and Leonardo Nicola Martin Iasi and Ekaterina Noskova and Jana Ob{\v{s}}teter and
+           Vitor Antonio Corr{\^{e}}a Pavinato and Alice Pearson and David Peede and Manolo F. Perez and
+           Murillo F. Rodrigues and Chris C. R. Smith and Jeffrey P. Spence and Anastasia Teterina and
+           Silas Tittes and Per Unneberg and Juan Manuel Vazquez and Ryan K. Waples and Anthony Wilder Wohns and
+           Yan Wong and Franz Baumdicker and Reed A. Cartwright and Gregor Gorjanc and Ryan N. Gutenkunst and
+           Jerome Kelleher and Andrew D. Kern and Aaron P. Ragsdale and Peter L. Ralph and Daniel R. Schrider and
+           Ilan Gronau},
+	doi = {10.7554/elife.84874.2},
+	url = {https://doi.org/10.7554/elife.84874.2},
+	year = 2023,
+	month = {may},
+	publisher = {{eLife} Sciences Publications, Ltd},
     }
 
 
diff --git a/stdpopsim/catalog/PhoSin/__init__.py b/stdpopsim/catalog/PhoSin/__init__.py
new file mode 100644
index 000000000..a518c0337
--- /dev/null
+++ b/stdpopsim/catalog/PhoSin/__init__.py
@@ -0,0 +1,4 @@
+"""
+Catalog definitions for PhoSin (Ensembl ID='phocoena_sinus')
+"""
+from . import species  # noqa: F401
diff --git a/stdpopsim/catalog/PhoSin/genome_data.py b/stdpopsim/catalog/PhoSin/genome_data.py
new file mode 100644
index 000000000..decb72cfb
--- /dev/null
+++ b/stdpopsim/catalog/PhoSin/genome_data.py
@@ -0,0 +1,29 @@
+# File autogenerated from Ensembl REST API. Do not edit.
+data = {
+    "assembly_accession": "GCA_008692025.1",
+    "assembly_name": "mPhoSin1.pri",
+    "chromosomes": {
+        "1": {"length": 185845356, "synonyms": []},
+        "2": {"length": 178563925, "synonyms": []},
+        "3": {"length": 174291665, "synonyms": []},
+        "4": {"length": 146127150, "synonyms": []},
+        "5": {"length": 139762554, "synonyms": []},
+        "6": {"length": 115952311, "synonyms": []},
+        "7": {"length": 115469292, "synonyms": []},
+        "8": {"length": 110408561, "synonyms": []},
+        "9": {"length": 106736052, "synonyms": []},
+        "10": {"length": 102828027, "synonyms": []},
+        "11": {"length": 104118372, "synonyms": []},
+        "12": {"length": 90399378, "synonyms": []},
+        "13": {"length": 90400161, "synonyms": []},
+        "14": {"length": 89762611, "synonyms": []},
+        "15": {"length": 88287594, "synonyms": []},
+        "16": {"length": 85252673, "synonyms": []},
+        "17": {"length": 79603858, "synonyms": []},
+        "18": {"length": 79885847, "synonyms": []},
+        "19": {"length": 59501331, "synonyms": []},
+        "20": {"length": 58878645, "synonyms": []},
+        "21": {"length": 35802323, "synonyms": []},
+        "X": {"length": 131664873, "synonyms": []},
+    },
+}
diff --git a/stdpopsim/catalog/PhoSin/species.py b/stdpopsim/catalog/PhoSin/species.py
new file mode 100644
index 000000000..6d3bd3fde
--- /dev/null
+++ b/stdpopsim/catalog/PhoSin/species.py
@@ -0,0 +1,173 @@
+import stdpopsim
+
+from . import genome_data
+
+_species_ploidy = 2
+_ploidy = {
+    "1": _species_ploidy,
+    "2": _species_ploidy,
+    "3": _species_ploidy,
+    "4": _species_ploidy,
+    "5": _species_ploidy,
+    "6": _species_ploidy,
+    "7": _species_ploidy,
+    "8": _species_ploidy,
+    "9": _species_ploidy,
+    "10": _species_ploidy,
+    "11": _species_ploidy,
+    "12": _species_ploidy,
+    "13": _species_ploidy,
+    "14": _species_ploidy,
+    "15": _species_ploidy,
+    "16": _species_ploidy,
+    "17": _species_ploidy,
+    "18": _species_ploidy,
+    "19": _species_ploidy,
+    "20": _species_ploidy,
+    "21": _species_ploidy,
+    "X": _species_ploidy,
+}
+
+##########################
+# There is no direct estimate of recombination rate
+# for Vaquita or any closely related species.
+# Other studies that simulate data for Vaquita assume
+# a default mutation rate of 1e-8, which is what we
+# use here.
+# For example, Morin et al. (2021) assume a recombination
+# rate of 1e-8 when running PSMC, and Robinson et al. (2022)
+# use a recombination rate of 1e-8 when simulating data using
+# SLiM.
+##########################
+_overall_recombination_rate = 1e-8
+_recombination_rate = {
+    "1": _overall_recombination_rate,
+    "2": _overall_recombination_rate,
+    "3": _overall_recombination_rate,
+    "4": _overall_recombination_rate,
+    "5": _overall_recombination_rate,
+    "6": _overall_recombination_rate,
+    "7": _overall_recombination_rate,
+    "8": _overall_recombination_rate,
+    "9": _overall_recombination_rate,
+    "10": _overall_recombination_rate,
+    "11": _overall_recombination_rate,
+    "12": _overall_recombination_rate,
+    "13": _overall_recombination_rate,
+    "14": _overall_recombination_rate,
+    "15": _overall_recombination_rate,
+    "16": _overall_recombination_rate,
+    "17": _overall_recombination_rate,
+    "18": _overall_recombination_rate,
+    "19": _overall_recombination_rate,
+    "20": _overall_recombination_rate,
+    "21": _overall_recombination_rate,
+    "X": _overall_recombination_rate,
+}
+
+##########################
+# Genome-wide average mutation rate is estimated by Robinson et al. (2022)
+# by using divergence from two closely related species: harbor porpoise
+# (Phocoena phocoena) and Indo-Pacific finless porpoise (Neophocaena phocaenoides).
+# Using these two species and a range of assumptions results in a range of
+# plausible estimates depicted in Fig S12 in the paper.
+# Eventually, the representative value used in most analyses in the paper was
+# 5.83e-9, which is on the high side of the range, to take into account previous
+# higher estimates of the mutation rate, cited by Morin et al. (2021).
+##########################
+_overall_mutation_rate = 5.83e-9
+_mutation_rate = {
+    "1": _overall_mutation_rate,
+    "2": _overall_mutation_rate,
+    "3": _overall_mutation_rate,
+    "4": _overall_mutation_rate,
+    "5": _overall_mutation_rate,
+    "6": _overall_mutation_rate,
+    "7": _overall_mutation_rate,
+    "8": _overall_mutation_rate,
+    "9": _overall_mutation_rate,
+    "10": _overall_mutation_rate,
+    "11": _overall_mutation_rate,
+    "12": _overall_mutation_rate,
+    "13": _overall_mutation_rate,
+    "14": _overall_mutation_rate,
+    "15": _overall_mutation_rate,
+    "16": _overall_mutation_rate,
+    "17": _overall_mutation_rate,
+    "18": _overall_mutation_rate,
+    "19": _overall_mutation_rate,
+    "20": _overall_mutation_rate,
+    "21": _overall_mutation_rate,
+    "X": _overall_mutation_rate,
+}
+
+_genome = stdpopsim.Genome.from_data(
+    genome_data.data,
+    recombination_rate=_recombination_rate,
+    mutation_rate=_mutation_rate,
+    ploidy=_ploidy,
+    citations=[
+        stdpopsim.Citation(
+            author="Morin et al.",
+            year=2021,
+            doi="https://doi.org/10.1111/1755-0998.13284",
+            reasons={
+                stdpopsim.CiteReason.ASSEMBLY,
+                stdpopsim.CiteReason.REC_RATE,
+            },
+        ),
+        stdpopsim.Citation(
+            author="Robinson et al.",
+            year=2022,
+            doi="https://doi.org/10.1126/science.abm1742",
+            reasons={stdpopsim.CiteReason.MUT_RATE},
+        ),
+    ],
+)
+
+_species = stdpopsim.Species(
+    id="PhoSin",
+    ensembl_id="phocoena_sinus",
+    name="Phocoena sinus",
+    common_name="Vaquita",
+    genome=_genome,
+    ##########################
+    # Robinson et al. (2022) and Morin et al. (2021) assumed an average
+    # generation time of 11.9 years in their analyses/calibrations.
+    # They cite Taylor et al. (2007) as a primary source.
+    # Table 1 in that paper specifies estimates of generation
+    # times for 58 cetacean species. It seems like the value selected was
+    # the one inferred for the sister species of harbor porpoise (Phocoena
+    # phocoena). The value estimated for Vaquita is slightly smaller (11.4),
+    # but maybe it reflects a recent reducion in generation time (???)
+    ##########################
+    generation_time=11.9,
+    ##########################
+    # Demographic models inferred for Vaquita by Morin et al. (2021) and
+    # Robinson et al. (2022) typically estimate small Ne of 2,000 - 5,000
+    # in the past 500,000 years, likely followed by some sharp bottleneck.
+    # We selected as a representative Ne a value of 3500, which is also
+    # consistent with the average Ne inferred in the 2-epoch model of
+    # Robinson et al. (2022)
+    ##########################
+    population_size=3500,
+    citations=[
+        stdpopsim.Citation(
+            author="Robinson et al.",
+            year=2022,
+            doi="https://doi.org/10.1126/science.abm1742",
+            reasons={
+                stdpopsim.CiteReason.POP_SIZE,
+                stdpopsim.CiteReason.GEN_TIME,
+            },
+        ),
+        stdpopsim.Citation(
+            author="Taylor et al.",
+            year=2007,
+            doi="https://aquadocs.org/handle/1834/41281",
+            reasons={stdpopsim.CiteReason.GEN_TIME},
+        ),
+    ],
+)
+
+stdpopsim.register_species(_species)


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 770b06160..bef554f7c 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -48,7 +48,7 @@ jobs:
         uses: actions/cache@v3
         env:
           # Increase this to reset the cache if the key hasn't changed.
-          CACHE_NUM: 5
+          CACHE_NUM: 6
         with:
           path: |
             ${{ steps.find-conda.outputs.CONDA }}/envs/${{ env.CONDA_ENV_NAME }}
diff --git a/tests/test_PhoSin.py b/tests/test_PhoSin.py
new file mode 100644
index 000000000..be93fe6ff
--- /dev/null
+++ b/tests/test_PhoSin.py
@@ -0,0 +1,99 @@
+import pytest
+
+import stdpopsim
+from tests import test_species
+
+
+class TestSpeciesData(test_species.SpeciesTestBase):
+
+    species = stdpopsim.get_species("PhoSin")
+
+    def test_ensembl_id(self):
+        assert self.species.ensembl_id == "phocoena_sinus"
+
+    def test_name(self):
+        assert self.species.name == "Phocoena sinus"
+
+    def test_common_name(self):
+        assert self.species.common_name == "Vaquita"
+
+    # QC Tests. These tests are performed by another contributor
+    # independently referring to the citations provided in the
+    # species definition, filling in the appropriate values
+    # and deleting the pytest "skip" annotations.
+    @pytest.mark.skip("Population size QC not done yet")
+    def test_qc_population_size(self):
+        assert self.species.population_size == -1
+
+    @pytest.mark.skip("Generation time QC not done yet")
+    def test_qc_generation_time(self):
+        assert self.species.generation_time == -1
+
+
+class TestGenomeData(test_species.GenomeTestBase):
+
+    genome = stdpopsim.get_species("PhoSin").genome
+
+    @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,
+            "X": -1,
+        }.items(),
+    )
+    def test_recombination_rate(self, name, rate):
+        assert rate == pytest.approx(
+            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,
+            "X": -1,
+        }.items(),
+    )
+    def test_mutation_rate(self, name, rate):
+        assert rate == pytest.approx(self.genome.get_chromosome(name).mutation_rate)