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

Issue:
QC for HomSap
The **HomSap** species definition has not yet been QCed!

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).

- [x] 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.
- [x] 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.
- [x] 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?
- [x] Population size.
- [x] 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:
- [x] Genome reference assembly.
- [x] Mutation rate.
- [x] Recombination rate.
- [x] Recombination map(s) (if relevant).
- [x] Population size.
- [x] Generation time.



Discussion / hints:
Over in #914 I'm finding the recombination rates *don't* match, eg:
```
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[13] - AssertionError: assert 1.0914860554958317e-08 == 1.3082718764398148e-08 ± 1.3e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[14] - AssertionError: assert 1.119730771394731e-08 == 1.3711030325976108e-08 ± 1.4e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[15] - AssertionError: assert 1.3835785893339787e-08 == 1.7197030015823156e-08 ± 1.7e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[17] - AssertionError: assert 1.582489036239487e-08 == 1.5435336515997285e-08 ± 1.5e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[18] - AssertionError: assert 1.5075956950023575e-08 == 1.4647379553777267e-08 ± 1.5e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[21] - AssertionError: assert 1.3045214034879191e-08 == 1.6862572870496712e-08 ± 1.7e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[X] - AssertionError: assert 1.164662223273842e-08 == 1.1793445889048436e-08 ± 1.2e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[20] - AssertionError: assert 1.7178269031631664e-08 == 1.6816238743401626e-08 ± 1.7e-10
FAILED tests/test_HomSap.py::TestGenome::test_recombination_rates[22] - AssertionError: assert 1.4445022767788226e-08 == 2.1066745562340675e-08 ± 2.1e-10
```
... so perhaps that first box should be un-checked.
Done! But I think these tests are going to be problematic because of the current vagueries associated with calculating rate averages across maps (we've changed how this could work several times in the past 12 months).
This one is up for grabs - and, it's almost done!
This one is not so easy, as it requires dealing with recombination map liftovers.
Ah. Darn it. Well, we've got docs for this: https://popsim-consortium.github.io/stdpopsim-docs/latest/development.html#lifting-over-a-genetic-map
I've worked with recombination map liftovers before. Happy to work on this. 
Did you get this lifted over, @saurabhbelsare?
I've shared the lifted over files with @andrewkern to upload to AWS, and was going to submit a pull request when he gave me the new updated URL. Sorry I forgot to follow up with him about this. 

Fix patch:
diff --git a/stdpopsim/catalog/HomSap/species.py b/stdpopsim/catalog/HomSap/species.py
index db831416a..4ce982874 100644
--- a/stdpopsim/catalog/HomSap/species.py
+++ b/stdpopsim/catalog/HomSap/species.py
@@ -2,31 +2,33 @@
 
 from . import genome_data
 
-# mean_recombination_rate was computed across all windows of the HapMapII_GRCh37 map
+# Mean recombination rate was computed from the "HapMapII_GRCh38" map
+# by dividing total map length by total sequence length over non-missing
+# windows. (See: `msprime.RateMap.mean_rate`)
 _recombination_rate_data = {
-    "1": 1.1485597641285933e-08,
-    "2": 1.1054289277533446e-08,
-    "3": 1.1279585624662551e-08,
-    "4": 1.1231162636001008e-08,
-    "5": 1.1280936570022824e-08,
-    "6": 1.1222852661225285e-08,
-    "7": 1.1764614397655721e-08,
-    "8": 1.1478465778920576e-08,
-    "9": 1.1780701596308656e-08,
-    "10": 1.3365134257075317e-08,
-    "11": 1.1719334320833283e-08,
-    "12": 1.305017186986983e-08,
-    "13": 1.0914860554958317e-08,
-    "14": 1.119730771394731e-08,
-    "15": 1.3835785893339787e-08,
-    "16": 1.4834607113882717e-08,
-    "17": 1.582489036239487e-08,
-    "18": 1.5075956950023575e-08,
-    "19": 1.8220141872466202e-08,
-    "20": 1.7178269031631664e-08,
-    "21": 1.3045214034879191e-08,
-    "22": 1.4445022767788226e-08,
-    "X": 1.164662223273842e-08,
+    "1": 1.1523470111585671e-08,
+    "2": 1.1042947599187705e-08,
+    "3": 1.1258496243179589e-08,
+    "4": 1.148200397734684e-08,
+    "5": 1.1244324085272748e-08,
+    "6": 1.1265873463291707e-08,
+    "7": 1.177126324439884e-08,
+    "8": 1.1604882118427952e-08,
+    "9": 1.219866363355389e-08,
+    "10": 1.3333710355319759e-08,
+    "11": 1.1721276661036567e-08,
+    "12": 1.3098075878935492e-08,
+    "13": 1.3061005794618377e-08,
+    "14": 1.3629801391728677e-08,
+    "15": 1.738760443308643e-08,
+    "16": 1.4831470988338458e-08,
+    "17": 1.55382549450431e-08,
+    "18": 1.4645528336185072e-08,
+    "19": 1.8384797447801312e-08,
+    "20": 1.6788552013988714e-08,
+    "21": 1.724434082825555e-08,
+    "22": 2.1057233894035443e-08,
+    "X": 1.1848323598825364e-08,
     "Y": 0.0,
     "MT": 0.0,
 }


Test patch:
diff --git a/tests/test_HomSap.py b/tests/test_HomSap.py
index 1b85a9f35..85c9e369d 100644
--- a/tests/test_HomSap.py
+++ b/tests/test_HomSap.py
@@ -18,34 +18,22 @@ def test_basic_attributes(self):
 
 class TestGenome(test_species.GenomeTestBase):
 
-    genome = stdpopsim.get_species("HomSap").genome
+    species = stdpopsim.get_species("HomSap")
+    genome = species.genome
 
     def test_basic_attributes(self):
         assert len(self.genome.chromosomes) == 25
 
-    @pytest.mark.skip("Failing - due to mismatching assembly?")
     @pytest.mark.parametrize("chr_id", [chrom.id for chrom in genome.chromosomes])
     def test_recombination_rates(self, chr_id):
-        # We should recast this test and just hard code in the values.
-        # Tests should be *obvious* not clever.
-
-        # recompute recombination rates from HapMapII_GRCh37 map then
-        # compare the results to the current recombination rates for each chromosome
-        genetic_map = "HapMapII_GRCh37"
-        species = stdpopsim.get_species("HomSap")
-        chrom = species.genome.get_chromosome(chr_id)
+        genetic_map = "HapMapII_GRCh38"
+        chrom = self.genome.get_chromosome(chr_id)
         if chr_id in ["X", "Y", "MT"]:
             with pytest.warns(stdpopsim.NonAutosomalWarning):
-                contig = species.get_contig(chr_id, genetic_map=genetic_map)
-        elif chr_id in ["3", "5", "7", "11", "16", "17", "18", "20"]:
-            contig = species.get_contig(chr_id, genetic_map=genetic_map)
+                contig = self.species.get_contig(chr_id, genetic_map=genetic_map)
         else:
-            # The rest of the chromosomes are currently emitting a warning about
-            # the mismatch in chromosome lengths because of the fact that we're
-            # on 37 for the map. This should be resolved when we start using the
-            # lifted over map.
-            with pytest.warns(UserWarning, match="longer than chromosome length"):
-                contig = species.get_contig(chr_id, genetic_map=genetic_map)
+            contig = self.species.get_contig(chr_id, genetic_map=genetic_map)
+        assert chrom.length == contig.recombination_map.sequence_length
         assert chrom.recombination_rate == pytest.approx(
             contig.recombination_map.mean_rate, rel=1e-6
         )
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index e7ceebd39..c879bb19b 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -942,7 +942,7 @@ def get_example_dfes(self):
     def test_slim_simulates_dfes(self):
         contig = get_test_contig()
         # make theta = 40
-        contig.mutation_rate = 10 / (1000 * contig.recombination_map.sequence_length)
+        contig.mutation_rate = 20 / (1000 * contig.recombination_map.sequence_length)
         print(contig.mutation_rate)
         dfes = [
             stdpopsim.DFE(