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

Issue:
Add Chlamydomonas reinhardtii to catalogue
- **Name:** _Chlamydomonas reinhardtii_
- **Common name:** n/a

**Description:** unicellular isogamous alga with a haploid 110 Mb genome - used in studies of cell motility, but also more recently as a model for evolutionary genetics (especially in light of its facultatively sexual life cycle) 

---

*Note: Maybe you don't know a lot of the things below - but, that's OK!
The criteria are that (a) the choice is reasonable, and (b) the rationale is clear.
What we want is a "reasonable estimate" (or, "best guess") - if you were to use a number in a paper,
what would you use and how would you justify it to your collaborators?
For instance, for "generation time", best would be to cite a paper that actually measures
or estimates generation time. But most species don't have this;
next best would be to give a number used in the literature, and provide a citation that used it
(and hopefully the publication gives a justification also).
Lacking this, you might provide a number from a related species (and a citation for this).*

*Here is the checklist of things that we need to add a species to the catalogue.
Each thing should be provided, with a justification (maybe short) and a citation.*

**Demographic information:**

- [x] generation time in years
  - ~2.4 mitotic doublings per day -> 365 * 2.4 = 876 doublings per year (Vitova et al 2011, [link][vitova2011])
  - 1/876 = 0.0011 years/generation
  - I hope I did this right - please correct me if I'm wrong! 
- [x] "default" population size
  - 1.4 * 10^7 (Ness et al 2016, [link](https://academic.oup.com/mbe/article/33/3/800/2579627?login=true))

[vitova2011]: https://link.springer.com/article/10.1007/s00425-011-1427-7

**Chromosome structure:**

*Note:* The assembly should be *chromosome-level*, i.e., not composed of thousands of scaffolds.

- [x] ensembl assembly ID *or*, if the assembly is not in Ensembl, a list of chromosomes with *name* and *length* (in bp)
  - GCA_000002595.3 ([Ensembl (Plants) link][ensembl], [ENA link][ena])

[ensembl]: http://plants.ensembl.org/Chlamydomonas_reinhardtii/Info/Annotation/#assembly
[ena]: https://www.ebi.ac.uk/ena/browser/view/GCA_000002595.3

**Recombination rates:**

- [ ] genetic map of recombination rates (as a hapmap or csv file)
- [x] genome-wide mean recombination rate (which will be used as the default)
  - 0.015 cM/Mb (Hasan and Ness 2020, [link][ld-paper])
  - This is = rho / 2Ne, and is two orders of magnitude lower than [lab estimates][kathir2003] (9.15 cM/Mb)! But if simulating wild populations of C. reinhardtii, in the absence of a means of modifying the rate of sex in natural populations, this value already takes the infrequency of sex into account
  - If a means of modifying the rate of sex is implemented into `stdpopsim`, then the 9.15 cM/Mb value is appropriate to use when simulating natural populations, provided the rate of sex is set to 1/840 (see [Hasan and Ness 2020][ld-paper])
  - But it is important to note that this is NOT the expected recombination rate in a single cross of C. reinhardtii, and is in fact substantially lower! 

[kathir2003]: https://ec.asm.org/content/2/2/362.short
[ld-paper]: https://academic.oup.com/gbe/article/12/4/370/5809146?login=true

**Mutation rate:**

- [x] genome-wide mean mutation rate
  - μ = 9.63 × 10^−10 (Ness et al 2015, [link][ness-2015])

[ness-2015]: https://genome.cshlp.org/content/25/11/1739.short


**Demographic model:**

- [ ] as a list of population sizes, growth rates, migration rates, etcetera. **(optional)**
- [ ] citation


**Other information:**

These are things we don't currently use, but will want to use in the future:

- [x] sex determination system, and which chromosomes are the sex chromosome(s)
  - C. reinhardtii is an isogamous, heterothallic species with mating types determined by a ~500 kb locus on chromosome 6
  - individuals are either MT+ or MT- based on which allele is present at said locus
  - see [Umen 2011][umen-2011] for more details
- [x] ploidy
  - haploid

[umen-2011]: www.sciencedirect.com/science/article/pii/S1369527411001755


Discussion / hints:

Fix patch:
diff --git a/stdpopsim/catalog/ChlRei/__init__.py b/stdpopsim/catalog/ChlRei/__init__.py
new file mode 100644
index 000000000..a6126ecfb
--- /dev/null
+++ b/stdpopsim/catalog/ChlRei/__init__.py
@@ -0,0 +1,4 @@
+"""
+Catalog definitions for ChlRei (Ensembl ID='chlamydomonas_reinhardtii')
+"""
+from . import species  # noqa: F401
diff --git a/stdpopsim/catalog/ChlRei/genome_data.py b/stdpopsim/catalog/ChlRei/genome_data.py
new file mode 100644
index 000000000..2ac80adef
--- /dev/null
+++ b/stdpopsim/catalog/ChlRei/genome_data.py
@@ -0,0 +1,24 @@
+# File autogenerated from Ensembl REST API. Do not edit.
+data = {
+    "assembly_accession": "GCA_000002595.3",
+    "assembly_name": "Chlamydomonas_reinhardtii_v5.5",
+    "chromosomes": {
+        "1": {"length": 8033585, "synonyms": []},
+        "2": {"length": 9223677, "synonyms": []},
+        "3": {"length": 9219486, "synonyms": []},
+        "4": {"length": 4091191, "synonyms": []},
+        "5": {"length": 3500558, "synonyms": []},
+        "6": {"length": 9023763, "synonyms": []},
+        "7": {"length": 6421821, "synonyms": []},
+        "8": {"length": 5033832, "synonyms": []},
+        "9": {"length": 7956127, "synonyms": []},
+        "10": {"length": 6576019, "synonyms": []},
+        "11": {"length": 3826814, "synonyms": []},
+        "12": {"length": 9730733, "synonyms": []},
+        "13": {"length": 5206065, "synonyms": []},
+        "14": {"length": 4157777, "synonyms": []},
+        "15": {"length": 1922860, "synonyms": []},
+        "16": {"length": 7783580, "synonyms": []},
+        "17": {"length": 7188315, "synonyms": []},
+    },
+}
diff --git a/stdpopsim/catalog/ChlRei/species.py b/stdpopsim/catalog/ChlRei/species.py
new file mode 100644
index 000000000..5cdd435fd
--- /dev/null
+++ b/stdpopsim/catalog/ChlRei/species.py
@@ -0,0 +1,105 @@
+import stdpopsim
+
+from . import genome_data
+
+# IMPORTANT: these are recombination rate values with the rate of sex
+# incorporated - e.g. they are lower than the physical recombination rates on account
+# of being the 'realized' recombination rates in nature, since C. reinhardtii only
+# undergoes sex/recombination once every ~840 generations (Hasan and Ness 2020).
+# when simulating natural populations with these values, the infrequency
+# of sex should be implicitly taken into account. otherwise, to get estimates
+# of physical rates, divide all values by 0.001194. see Discussion of Hasan
+# and Ness 2020 for more details
+_recombination_rate = {
+    "1": 1.21e-10,
+    "2": 1.49e-10,
+    "3": 1.52e-10,
+    "4": 1.47e-10,
+    "5": 1.70e-10,
+    "6": 1.17e-10,
+    "7": 8.66e-11,
+    "8": 1.39e-10,
+    "9": 1.12e-10,
+    "10": 1.97e-10,
+    "11": 1.63e-10,
+    "12": 9.15e-11,
+    "13": 1.43e-10,
+    "14": 1.9e-10,
+    "15": 3.93e-10,
+    "16": 1.71e-10,
+    "17": 1.83e-10,
+}
+
+_overall_rate = 9.63 * 1e-10  # SNM rate
+# there is also an indel mutation rate in Ness 2015 if that is eventually supported
+_mutation_rate = {
+    "1": _overall_rate,
+    "2": _overall_rate,
+    "3": _overall_rate,
+    "4": _overall_rate,
+    "5": _overall_rate,
+    "6": _overall_rate,
+    "7": _overall_rate,
+    "8": _overall_rate,
+    "9": _overall_rate,
+    "10": _overall_rate,
+    "11": _overall_rate,
+    "12": _overall_rate,
+    "13": _overall_rate,
+    "14": _overall_rate,
+    "15": _overall_rate,
+    "16": _overall_rate,
+    "17": _overall_rate,
+}
+
+_genome = stdpopsim.Genome.from_data(
+    genome_data.data,
+    recombination_rate=_recombination_rate,
+    mutation_rate=_mutation_rate,
+    citations=[
+        stdpopsim.Citation(
+            author="Merchant et al",
+            year=2007,
+            doi="https://doi.org/10.1126/science.1143609",
+            reasons={stdpopsim.CiteReason.ASSEMBLY},  # v5 - v6 assembly still en route!
+        ),
+        stdpopsim.Citation(
+            author="Hasan and Ness",
+            year=2020,
+            doi="https://doi.org/10.1093/gbe/evaa057",
+            reasons={stdpopsim.CiteReason.REC_RATE},
+        ),
+        stdpopsim.Citation(
+            author="Ness et al",
+            year=2015,
+            doi="https://doi.org/10.1101/gr.191494.115",
+            reasons={stdpopsim.CiteReason.MUT_RATE},
+        ),
+    ],
+)
+
+_species = stdpopsim.Species(
+    id="ChlRei",
+    ensembl_id="chlamydomonas_reinhardtii",
+    name="Chlamydomonas reinhardtii",
+    common_name="Chlamydomonas reinhardtii",
+    genome=_genome,
+    generation_time=1 / 876,
+    population_size=1.4 * 1e-7,
+    citations=[
+        stdpopsim.Citation(
+            author="Ness et al.",
+            year=2016,
+            doi="https://doi.org/10.1093/molbev/msv272",
+            reasons={stdpopsim.CiteReason.POP_SIZE},  # Quebec population
+        ),
+        stdpopsim.Citation(
+            author="Vitova et al.",
+            year=2011,
+            doi="https://doi.org/10.1007/s00425-011-1427-7",
+            reasons={stdpopsim.CiteReason.GEN_TIME},
+        ),
+    ],
+)
+
+stdpopsim.register_species(_species)


Test patch:
diff --git a/tests/test_ChlRei.py b/tests/test_ChlRei.py
new file mode 100644
index 000000000..5bf36aa98
--- /dev/null
+++ b/tests/test_ChlRei.py
@@ -0,0 +1,87 @@
+import pytest
+
+import stdpopsim
+from tests import test_species
+
+
+class TestSpeciesData(test_species.SpeciesTestBase):
+
+    species = stdpopsim.get_species("ChlRei")
+
+    def test_ensembl_id(self):
+        assert self.species.ensembl_id == "chlamydomonas_reinhardtii"
+
+    def test_name(self):
+        assert self.species.name == "Chlamydomonas reinhardtii"
+
+    def test_common_name(self):
+        assert self.species.common_name == "Chlamydomonas reinhardtii"
+
+    # 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("ChlRei").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,
+        }.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,
+        }.items(),
+    )
+    def test_mutation_rate(self, name, rate):
+        assert pytest.approx(rate, self.genome.get_chromosome(name).mutation_rate)