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

Issue:
QC for Oryza sativa
**PR for new species:** Oryza sativa #1453 

If you volunteer to QC this species, please use the checklist below.
While this list is intended to be comprehensive, it may not be exhaustive.
Where relevant, the QC reviewer should identify that parameter values match
those given in the linked citation(s).

The QC reviewer should start a pull request that fills out the test stubs
with independently obtained values. The reviewer may look at the python code
for rationale provided in comments, but should ignore the actual code
as much as possible - comments in the code should give enough information
that it's obvious how to get the correct value from the provided references.
(In particular, we shouldn't copy-paste the value from the code into the test!)

- [ ] Recombination rate.
  - This might be genome-wide, or per-chromosome. Both are fine.
  - Check there's a comment describing where it came from, and/or how calculated.
  - From a publication? Check the value(s) match the publication.
  - Calculated somehow? Average over a recombination map? Redo the calculation.
- [ ] Mutation rate.
  - This might be genome-wide, or per-chromosome. Both are fine.
  - Check there's a comment describing where it came from, and/or how calculated.
  - From a publication? Check the value(s) match the publication.
  - Calculated somehow? Redo the calculation.
- [ ] Recombination map (if present).
  - Does it match the assembly? Liftover is fine, if clearly stated.
  - Is the description/long_description a good summary of how the map was created?
- [ ] Population size.
- [ ] Generation time.

For each citation, check:
- Doi link.
- Is publication a preprint? Is there a peer-reviewed publication instead?
- Is the year correct.
- Is the author correct (spelling, accents/ligatures/etc.)

Citations are required for:
- [ ] Genome reference assembly.
- [ ] Mutation rate.
- [ ] Recombination rate.
- [ ] Recombination map(s) (if relevant).
- [ ] Population size.
- [ ] Generation time.

The final PR should:
- [ ] fill out the test stubs
- [ ] delete the `pytest.mark.skip` lines that make the tests not run
- [ ] make sure they pass, talking to the original author to figure out discrepancies



Discussion / hints:
Thanks for volunteering to QC, @minesrebollo!
Hi! 

This is my progress so far on checking this:

- [x] Recombination rate
- [x] Genome assembly
- [x] Generation time

 I don't find the provided mutation rate number in the cited reference. I had previously used 3.2e-9 from [Yang et al. 2015](https://www.nature.com/articles/nature14649), what do you think about using this one instead @ornobalam? Or maybe you can point out where the provided number is on Gaut et al. 1996. 

Citations are checked for:
- [x] Genome reference assembly
- [x] Recombination rate
- [x] Generation time

I will check the population size and demographic model later this week, sorry for the delay on doing this!
Hi! So, I picked up this mutation rate and citation from a few different papers (including [Salse et al. 2008](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2254919/), [Liu et al. 2017](https://academic.oup.com/mbe/article/34/4/908/2897207), and [Wang et al. 2016](https://www.sciencedirect.com/science/article/pii/S1674205216300557)). My best guess for how people calculated it comes from this passage in Gaut et al. 1996:
"The grass Adhl locus is estimated to evolve at a rate of 7.00 x 10-9 (±1.16 x 10-9) substitution per synonymous site per year, and the Adh2 locus is estimated to evolve at a rate of 5.99 x 10-9 (±1.13 x 10-9) substitution per synonymous site per year (Fig. 3A)."

They probably just calculated the mean of the rates for the two loci. Should've added all this as a comment!

I hadn't seen Yang et al. 2015, but I think it's clearly more robust. And using the 3.2e-9 rate for homozygotes makes sense as you see little heterozygosity in rice.

Thanks for going through these! 
Hi @ornobalam, your guess on how those authors got the recombination rate from Gaut et al. 1996 sounds good but if you agree that Yang et al. 2015 is more robust I think that is better too. Will you make that change in the commit or will it be part of my review? I am not sure about this process... 

For the population size, I see that you used the ratios of the effective population size of the current population to ancestral for the model with migration + bottleneck (v2 in table 2) but I wonder about the ancestral (rufipogon) population size estimate. They provide three different estimates, one from indica analyses, Nrufi = 57,471, other from tropical japonica, Nrufi = 118,110 and an independent one using the estimated scaled population silent mutation rates, where Nrufi = 135,500 but none is the one you used (150,000). Could you please point out what might I be missing? :)

Thanks!
I now see that you specified this in `N_rufi = 12000 / (0.04 * 2)`, I wonder why did they not present this estimate in the paper too.
Migration rates look good to me.

@nspope, after clarifying where the change in the mutation rate will be done and agreeing that the N_rufi used is the right one, do I need to start a pull request or that is all?

Hi! @minesrebollo, I am not sure about the process either, but I can change the mutation rate.

I agree it's not clear why they didn't present the 150,000 estimate in the paper. Especially because the migration model produced better fits for the different SFSs than the individual bottleneck models, which is why I used this estimate. My reason for not using the 135,000 estimate was that they calculated it from theta using a mutation rate of 1e-8.

There are of course other papers that provide Ne estimates for rufipogon. For instance, [Molina et al. 2011](https://www.pnas.org/doi/10.1073/pnas.1104686108) provides two estimates (~450,000 and ~260,000, both found in the supplementary information pdf). These are based on calculating theta from the number of segregating sites and then using different mutation rates to calculate Ne from theta=4Neu. It's not clear to me which estimate may be better across studies, but I do think the 150,000 estimate makes the most sense to use from Caicedo et al.

(I have also been wondering about something else; the population size to be provided is described in the documentation like this: "the effective population size should represent the historical average effective population size". I interpreted this to mean ancestral effective population. Is that fair, @nspope ? The effective population sizes of the two rice subspecies japonica and indica (which are quite diverged but with gene flow) are distinct and the only way to have a historical effective size is through considering wild ancestral rice.)

If my interpretation of the effective population size here is correct, I propose to use the 150,000 estimate from Caicedo et al.

Thanks!

> @nspope, after clarifying where the change in the mutation rate will be done and agreeing that the N_rufi used is the right one, do I need to start a pull request or that is all?

Yup -- please start a pull request that fills out the test stubs in `tests/test_OrySat.py`, and fixes parameters in the various `catalog/OrySat/*.py` files. Thank you!
> the population size to be provided is described in the documentation like this: "the effective population size should represent the historical average effective population size". I interpreted this to mean ancestral effective population. Is that fair, @nspope ?

Yeah, good question! Rather than ancestral population size, I think what is meant here is a Ne value that produces levels of nucleotide diversity similar to that found in real data. For example, if you have a species that has gone through a very recent population decline, using its "ancestral Ne" for a global population size might make simulations unrealistically diverse. Maybe @igronau can clarify?
I agree with @nspope - however, there's not a single correct answer; the idea is to pick a reasonable value and document how you made the choice. (So, if it's not easy to find a documented value for nucleotide diversity (plus mut rate) to transform into an Ne then taking some other number from the demographic model - either modern or ancestral size - could be fine also.)

Fix patch:
diff --git a/stdpopsim/catalog/OrySat/species.py b/stdpopsim/catalog/OrySat/species.py
index b98bfa73d..3ee54ff38 100644
--- a/stdpopsim/catalog/OrySat/species.py
+++ b/stdpopsim/catalog/OrySat/species.py
@@ -70,6 +70,7 @@
 )
 # Set some exceptions for non-recombining chrs.
 _recombination_rate_data["Mt"] = 0
+_recombination_rate_data["Pt"] = 0
 
 # Generic and chromosome-specific ploidy
 _species_ploidy = 2


Test patch:
diff --git a/tests/test_OrySat.py b/tests/test_OrySat.py
index 8da2782f3..7a82fd0d8 100644
--- a/tests/test_OrySat.py
+++ b/tests/test_OrySat.py
@@ -16,40 +16,33 @@ def test_name(self):
     def test_common_name(self):
         assert self.species.common_name == "Asian rice"
 
-    # 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
+        assert self.species.population_size == 46875
 
-    @pytest.mark.skip("Generation time QC not done yet")
     def test_qc_generation_time(self):
-        assert self.species.generation_time == -1
+        assert self.species.generation_time == 1
 
 
 class TestGenomeData(test_species.GenomeTestBase):
     genome = stdpopsim.get_species("OrySat").genome
 
-    @pytest.mark.skip("Recombination rate QC not done yet")
     @pytest.mark.parametrize(
         ["name", "rate"],
         {
-            "1": -1,
-            "Mt": -1,
-            "Pt": -1,
-            "2": -1,
-            "3": -1,
-            "4": -1,
-            "5": -1,
-            "6": -1,
-            "7": -1,
-            "8": -1,
-            "9": -1,
-            "10": -1,
-            "11": -1,
-            "12": -1,
+            "1": 8.97e-10,
+            "Mt": 0,
+            "Pt": 0,
+            "2": 8.97e-10,
+            "3": 8.97e-10,
+            "4": 8.97e-10,
+            "5": 8.97e-10,
+            "6": 8.97e-10,
+            "7": 8.97e-10,
+            "8": 8.97e-10,
+            "9": 8.97e-10,
+            "10": 8.97e-10,
+            "11": 8.97e-10,
+            "12": 8.97e-10,
         }.items(),
     )
     def test_recombination_rate(self, name, rate):
@@ -57,25 +50,31 @@ def test_recombination_rate(self, name, rate):
             self.genome.get_chromosome(name).recombination_rate
         )
 
-    @pytest.mark.skip("Mutation rate QC not done yet")
     @pytest.mark.parametrize(
         ["name", "rate"],
         {
-            "1": -1,
-            "Mt": -1,
-            "Pt": -1,
-            "2": -1,
-            "3": -1,
-            "4": -1,
-            "5": -1,
-            "6": -1,
-            "7": -1,
-            "8": -1,
-            "9": -1,
-            "10": -1,
-            "11": -1,
-            "12": -1,
+            "1": 3.2e-9,
+            "Mt": 3.2e-9,
+            "Pt": 3.2e-9,
+            "2": 3.2e-9,
+            "3": 3.2e-9,
+            "4": 3.2e-9,
+            "5": 3.2e-9,
+            "6": 3.2e-9,
+            "7": 3.2e-9,
+            "8": 3.2e-9,
+            "9": 3.2e-9,
+            "10": 3.2e-9,
+            "11": 3.2e-9,
+            "12": 3.2e-9,
         }.items(),
     )
     def test_mutation_rate(self, name, rate):
         assert rate == pytest.approx(self.genome.get_chromosome(name).mutation_rate)
+
+    @pytest.mark.parametrize("chrom", [chrom for chrom in genome.chromosomes])
+    def test_chromosome_ploidy(self, chrom):
+        if chrom.id in ["Mt", "Pt"]:
+            assert chrom.ploidy == 1
+        else:
+            assert chrom.ploidy == 2