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

Issue:
QC for BosTau
**PR for new species:** #600 #809 

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

QC for BosTau
**PR for new species:** #600 #809 

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:
I believe this hasn't been started - do you know otherwise, @gregorgorjanc?
I'll do this one!
Looking at [`tests/test_BosTau.py`](https://github.com/popsim-consortium/stdpopsim/blob/main/tests/test_BosTau.py) this *has* been started, but I couldn't find the issue or PR, and it looks like it's not finished - recombination and mutation rates [aren't filled in yet](https://github.com/popsim-consortium/stdpopsim/blob/main/tests/test_BosTau.py#L35). (I'm probably missing something?)
I'm having some difficulty finding the relevant information in the [Rosen et al. publication](https://academic.oup.com/gigascience/article/9/3/giaa021/5810242) that supports our BosTau genome assembly data -- in particular, the chromosome lengths and *any* MT information. Part of the supplementary materials (the genetic map referenced in S2) requires a publisher login -- perhaps it's in there? Could we get the original contributor in here for assistance, @petrelharp or @grahamgower?
@grahamgower mentions [here](https://github.com/popsim-consortium/stdpopsim/pull/600#issue-713518206) that he changed the species-wide population size (and maybe also generation time?) to reflect their higher ancestral value, rather than the small modern-day value. However, both the OG code and the beginning of the QC code use the smaller modern value of 90 -- has this been changed back in the interim, or is there something else I'm missing?
Re: generation time: at the moment, we cite [Macleod et al](https://doi.org/10.1093/molbev/mst125), but the authors themselves mention that they just chose this value to be in the same ballpark as previous attempts based on actual estimation and observation ([here]([10.2527/jas.2006-367](https://doi.org/10.2527/jas.2006-367)) and [here]([10.1186/1297-9686-35-1-43](https://dx.doi.org/10.1186%2F1297-9686-35-1-43))). I'm happy to just use the nice round heuristic value directly from Macleod et al, but it might be worth changing it later to once of these other references.
I don't recall the details here, sorry. Maybe @gregorgorjanc remembers? There are a bunch of related PRs and issues for BosTau, so maybe checking some of these will help? #558, #559, #579, #609
Thanks @grahamgower, looks like for [this issue at least](https://github.com/popsim-consortium/stdpopsim/issues/1239#issuecomment-1129495160), it's the citation that's wrong -- from #558 it sounds like that info's come from the Ensembl genome browser instead of the one currently lited. I'll update the citation as a separate PR
I believe this hasn't been started - do you know otherwise, @gregorgorjanc?
I'll do this one!
Looking at [`tests/test_BosTau.py`](https://github.com/popsim-consortium/stdpopsim/blob/main/tests/test_BosTau.py) this *has* been started, but I couldn't find the issue or PR, and it looks like it's not finished - recombination and mutation rates [aren't filled in yet](https://github.com/popsim-consortium/stdpopsim/blob/main/tests/test_BosTau.py#L35). (I'm probably missing something?)
I'm having some difficulty finding the relevant information in the [Rosen et al. publication](https://academic.oup.com/gigascience/article/9/3/giaa021/5810242) that supports our BosTau genome assembly data -- in particular, the chromosome lengths and *any* MT information. Part of the supplementary materials (the genetic map referenced in S2) requires a publisher login -- perhaps it's in there? Could we get the original contributor in here for assistance, @petrelharp or @grahamgower?
@grahamgower mentions [here](https://github.com/popsim-consortium/stdpopsim/pull/600#issue-713518206) that he changed the species-wide population size (and maybe also generation time?) to reflect their higher ancestral value, rather than the small modern-day value. However, both the OG code and the beginning of the QC code use the smaller modern value of 90 -- has this been changed back in the interim, or is there something else I'm missing?
Re: generation time: at the moment, we cite [Macleod et al](https://doi.org/10.1093/molbev/mst125), but the authors themselves mention that they just chose this value to be in the same ballpark as previous attempts based on actual estimation and observation ([here]([10.2527/jas.2006-367](https://doi.org/10.2527/jas.2006-367)) and [here]([10.1186/1297-9686-35-1-43](https://dx.doi.org/10.1186%2F1297-9686-35-1-43))). I'm happy to just use the nice round heuristic value directly from Macleod et al, but it might be worth changing it later to once of these other references.
I don't recall the details here, sorry. Maybe @gregorgorjanc remembers? There are a bunch of related PRs and issues for BosTau, so maybe checking some of these will help? #558, #559, #579, #609
Thanks @grahamgower, looks like for [this issue at least](https://github.com/popsim-consortium/stdpopsim/issues/1239#issuecomment-1129495160), it's the citation that's wrong -- from #558 it sounds like that info's come from the Ensembl genome browser instead of the one currently lited. I'll update the citation as a separate PR

Fix patch:
diff --git a/stdpopsim/catalog/BosTau/species.py b/stdpopsim/catalog/BosTau/species.py
index 948959bd9..086513807 100644
--- a/stdpopsim/catalog/BosTau/species.py
+++ b/stdpopsim/catalog/BosTau/species.py
@@ -11,6 +11,13 @@
     reasons={stdpopsim.CiteReason.ASSEMBLY},
 )
 
+_HoweEtAl = stdpopsim.Citation(
+    author="Howe et al.",
+    year=2020,
+    doi="https://doi.org/10.1093/nar/gkaa942",
+    reasons={stdpopsim.CiteReason.ASSEMBLY},
+)
+
 # Frequency of mosaicism points towards mutation-prone early cleavage
 # cell divisions in cattle.
 _HarlandEtAl = stdpopsim.Citation(
@@ -73,7 +80,7 @@
 _genome = stdpopsim.Genome(
     chromosomes=_chromosomes,
     citations=[
-        _RosenEtAl,
+        _HoweEtAl,
         _HarlandEtAl,
         _MaEtAl,
     ],
@@ -87,7 +94,7 @@
     common_name="Cattle",
     genome=_genome,
     generation_time=5,
-    population_size=90,  # most recent Ne in _MacLeodEtAl
+    population_size=62000,  # ancestral Ne in _MacLeodEtAl
     citations=[_MacLeodEtAl],
 )
 


Test patch:
diff --git a/tests/test_BosTau.py b/tests/test_BosTau.py
index c1e130a3a..c73e1e2ce 100644
--- a/tests/test_BosTau.py
+++ b/tests/test_BosTau.py
@@ -18,7 +18,7 @@ def test_common_name(self):
         assert self.species.common_name == "Cattle"
 
     def test_qc_population_size(self):
-        assert self.species.population_size == 90  # most recent Ne in _MacLeodEtAl
+        assert self.species.population_size == 62000  # ancestral Ne in _MacLeodEtAl
 
     def test_qc_generation_time(self):
         assert self.species.generation_time == 5
@@ -32,83 +32,128 @@ def test_basic_attributes(self):
         nchrom = 31  # 29 + X + MT
         assert len(self.genome.chromosomes) == nchrom
 
-    @pytest.mark.skip("Recombination rate QC not done yet")
+    @pytest.mark.parametrize(
+        ["name", "length"],
+        {
+            "1": 158534110,
+            "2": 136231102,
+            "3": 121005158,
+            "4": 120000601,
+            "5": 120089316,
+            "6": 117806340,
+            "7": 110682743,
+            "8": 113319770,
+            "9": 105454467,
+            "10": 103308737,
+            "11": 106982474,
+            "12": 87216183,
+            "13": 83472345,
+            "14": 82403003,
+            "15": 85007780,
+            "16": 81013979,
+            "17": 73167244,
+            "18": 65820629,
+            "19": 63449741,
+            "20": 71974595,
+            "21": 69862954,
+            "22": 60773035,
+            "23": 52498615,
+            "24": 62317253,
+            "25": 42350435,
+            "26": 51992305,
+            "27": 45612108,
+            "28": 45940150,
+            "29": 51098607,
+            "X": 139009144,
+            "MT": 16338,
+        }.items(),
+    )
+    def test_chromosome_lengths(self, name, length):
+        assert length == self.genome.get_chromosome(name).length
+
+    # Performs the calculation described in the species file comments
+    lengths = [c.length for c in genome.chromosomes]
+    total_length = sum(lengths) - genome.get_chromosome("MT").length
+    avg_crossovers = (25.5 + 23.2) / 2  # average of male and female numbers
+    avg_rate = avg_crossovers / total_length
+
     @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,
-            "22": -1,
-            "23": -1,
-            "24": -1,
-            "25": -1,
-            "26": -1,
-            "27": -1,
-            "28": -1,
-            "29": -1,
-            "X": -1,
-            "MT": -1,
+            "1": avg_rate,
+            "2": avg_rate,
+            "3": avg_rate,
+            "4": avg_rate,
+            "5": avg_rate,
+            "6": avg_rate,
+            "7": avg_rate,
+            "8": avg_rate,
+            "9": avg_rate,
+            "10": avg_rate,
+            "11": avg_rate,
+            "12": avg_rate,
+            "13": avg_rate,
+            "14": avg_rate,
+            "15": avg_rate,
+            "16": avg_rate,
+            "17": avg_rate,
+            "18": avg_rate,
+            "19": avg_rate,
+            "20": avg_rate,
+            "21": avg_rate,
+            "22": avg_rate,
+            "23": avg_rate,
+            "24": avg_rate,
+            "25": avg_rate,
+            "26": avg_rate,
+            "27": avg_rate,
+            "28": avg_rate,
+            "29": avg_rate,
+            "X": avg_rate,
+            "MT": 0.0,
         }.items(),
     )
     def test_recombination_rate(self, name, rate):
         assert rate == pytest.approx(
-            self.genome.get_chromosome(name).recombination_rate
+            self.genome.get_chromosome(name).recombination_rate, rel=0.001
         )
 
-    @pytest.mark.skip("Mutation rate QC not done yet")
+    _mutation_rate = 1.2e-8
+
     @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,
-            "22": -1,
-            "23": -1,
-            "24": -1,
-            "25": -1,
-            "26": -1,
-            "27": -1,
-            "28": -1,
-            "29": -1,
-            "X": -1,
-            "MT": 0,
+            "1": _mutation_rate,
+            "2": _mutation_rate,
+            "3": _mutation_rate,
+            "4": _mutation_rate,
+            "5": _mutation_rate,
+            "6": _mutation_rate,
+            "7": _mutation_rate,
+            "8": _mutation_rate,
+            "9": _mutation_rate,
+            "10": _mutation_rate,
+            "11": _mutation_rate,
+            "12": _mutation_rate,
+            "13": _mutation_rate,
+            "14": _mutation_rate,
+            "15": _mutation_rate,
+            "16": _mutation_rate,
+            "17": _mutation_rate,
+            "18": _mutation_rate,
+            "19": _mutation_rate,
+            "20": _mutation_rate,
+            "21": _mutation_rate,
+            "22": _mutation_rate,
+            "23": _mutation_rate,
+            "24": _mutation_rate,
+            "25": _mutation_rate,
+            "26": _mutation_rate,
+            "27": _mutation_rate,
+            "28": _mutation_rate,
+            "29": _mutation_rate,
+            "X": _mutation_rate,
+            "MT": _mutation_rate,
         }.items(),
     )
     def test_mutation_rate(self, name, rate):