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

Issue:
QC for AnoGam
**PR for new species:** #856 

refs used for the parameters/assembly follow:
https://doi.org/10.1186/gb-2007-8-1-r5
https://doi.org/10.4269/ajtmh.2006.75.901
https://doi.org/10.1038/nature24995

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).

- [ ] 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.



Discussion / hints:
Just updating here: the QC should include a PR that updates the ``tests/test_AnoGam.py`` to fill in the missing values and delete the "pytest.skip" lines. All these tests should pass before the species is QC'd.

The QC should be blind, so the QC person **should not** look at the catalog file.

Fix patch:
diff --git a/stdpopsim/catalog/AnoGam/species.py b/stdpopsim/catalog/AnoGam/species.py
index 39c003e31..821679d89 100644
--- a/stdpopsim/catalog/AnoGam/species.py
+++ b/stdpopsim/catalog/AnoGam/species.py
@@ -21,6 +21,14 @@
     reasons={stdpopsim.CiteReason.REC_RATE},
 )
 
+_ZhengEtAl = stdpopsim.Citation(
+    # there's no DOI?!?!
+    doi="https://www.genetics.org/content/143/2/941",
+    year=1996,
+    author="Zheng et al.",
+    reasons={stdpopsim.CiteReason.REC_RATE},
+)
+
 _SharakhovaEtAl = stdpopsim.Citation(
     doi="https://doi.org/10.1186/gb-2007-8-1-r5",
     year=2006,
@@ -28,17 +36,30 @@
     reasons={stdpopsim.CiteReason.ASSEMBLY},
 )
 
+_KeightleyEtAl = stdpopsim.Citation(
+    doi="https://doi.org/10.1101/gr.091231.109",
+    year=2009,
+    author="Keightley et al.",
+    reasons={stdpopsim.CiteReason.MUT_RATE},
+)
+
+# values are in Pombi et al, but except for the X,
+# they refer to Zheng et al for the data
+# and have no data for 2L so we set equal to 2R
+# the value for X is the average from Table 1 in Zheng
 _recombination_rate = {
-    "2L": 0,  # setting to zero because of inversion
+    "2L": 1.3e-8,
     "2R": 1.3e-8,
-    "3L": 1.6e-8,
-    "3R": 1.3e-8,
-    "X": 1e-8,
+    "3L": 1.3e-8,
+    "3R": 1.6e-8,
+    "X": 2.04e-8,
     "Mt": 0,
 }
 
 
-_overall_rate = 5.49e-9
+# the rate for DroMel as inferred by Keightley et al 2009,
+# which was used for the stairwayplot demographic inference
+_overall_rate = 3.5e-9
 _mutation_rate = {
     "2L": _overall_rate,
     "2R": _overall_rate,
@@ -56,6 +77,8 @@
         _SharakhovaEtAl,
         _Ag1000G,
         _PombiEtAl,
+        _KeightleyEtAl,
+        _ZhengEtAl,
     ],
 )
 stdpopsim.utils.append_common_synonyms(_genome)
@@ -67,7 +90,8 @@
     common_name="Anopheles gambiae",
     genome=_genome,
     generation_time=1 / 11,
-    population_size=6e6,  # Ghana population
+    # based on theta = 4 Ne u in Gabon population, rounded
+    population_size=1e6,
     citations=[_Ag1000G],
 )
 


Test patch:
diff --git a/tests/test_AnoGam.py b/tests/test_AnoGam.py
index d3356d189..72acf0bb1 100644
--- a/tests/test_AnoGam.py
+++ b/tests/test_AnoGam.py
@@ -17,37 +17,52 @@ def test_name(self):
     def test_common_name(self):
         assert self.species.common_name == "Anopheles gambiae"
 
-    # 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
+        # based on theta = 4 Ne u since diversity in most pops is 1.5%
+        assert self.species.population_size == pytest.approx(
+            0.015 / (4 * 3.5e-9), rel=0.1
+        )
 
-    @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 / 11
 
 
 class TestGenomeData(test_species.GenomeTestBase):
 
     genome = stdpopsim.get_species("AnoGam").genome
+    mut_rate = 3.5e-9
 
-    @pytest.mark.skip("Recombination rate QC not done yet")
+    # these values are given in the text of Pombi et al
+    # except no single value is given for the X; the number
+    # given is the average across chromosome from Pombi's Table 1
+    # No value is given for 2L, because of the inversion,
+    # so we set it equal to 2R.
     @pytest.mark.parametrize(
         ["name", "rate"],
-        {"2L": -1, "2R": -1, "3L": -1, "3R": -1, "X": -1, "Mt": -1}.items(),
+        {
+            "2L": 1.3e-8,
+            "2R": 1.3e-8,
+            "3L": 1.3e-8,
+            "3R": 1.6e-8,
+            "X": round(41.379 / (22.104923 - 1.825909) * 1e-8, 10),
+            "Mt": 0,
+        }.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"],
-        {"2L": -1, "2R": -1, "3L": -1, "3R": -1, "X": -1, "Mt": -1}.items(),
+        {
+            "2L": mut_rate,
+            "2R": mut_rate,
+            "3L": mut_rate,
+            "3R": mut_rate,
+            "X": mut_rate,
+            "Mt": mut_rate,
+        }.items(),
     )
     def test_mutation_rate(self, name, rate):
         assert rate == pytest.approx(self.genome.get_chromosome(name).mutation_rate)