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

Issue:
ZeroDivisionError during PonAbe liftover
Get necessary files.
```
$ wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver
$ wget https://hgdownload.soe.ucsc.edu/goldenPath/ponAbe3/liftOver/ponAbe3ToPonAbe2.over.chain.gz
$ wget https://hgdownload.soe.ucsc.edu/goldenPath/ponAbe2/liftOver/ponAbe2ToPonAbe3.over.chain.gz
```
Be sure to `chmod +x liftOver` and put it in your `$PATH`. Then try the `liftOver_catalog.py` script as outlined in the dev docs.

```
$ python maintenance/liftOver_catalog.py \
  --species PonAbe \
  --map NaterPA_PonAbe2 \
  --chainFile ponAbe2ToPonAbe3.over.chain.gz \
  --validationChain ponAbe3ToPonAbe2.over.chain.gz \
  --winLen 1000 \
  --gapThresh 1000000 \
  --useAdjacentAvg \
  --retainIntermediates

Downloading tarball...
Extracting tarball...

Mapping the rates to the new assembly...
Reading liftover chains
Mapping coordinates
Reading liftover chains
Reading liftover chains
Mapping coordinates
Mapping coordinates
Reading liftover chains
Reading liftover chains
Mapping coordinates
Mapping coordinates
Reading liftover chains
Mapping coordinates
Reading liftover chains
Reading liftover chains
Mapping coordinates
Mapping coordinates
Process Process-2:
Traceback (most recent call last):
  File "/usr/lib/python3.10/multiprocessing/process.py", line 315, in _bootstrap
    self.run()
  File "/usr/lib/python3.10/multiprocessing/process.py", line 108, in run
    self._target(*self._args, **self._kwargs)
  File "/home/grg/src/stdpopsim/maintenance/liftOver.py", line 561, in worker
    joinedRates, new_average = bed_to_hapmap(
  File "/home/grg/src/stdpopsim/maintenance/liftOver.py", line 143, in bed_to_hapmap
    rescaled_rate = rescale_rate(oldRateDict, ar)
  File "/home/grg/src/stdpopsim/maintenance/liftOver.py", line 124, in rescale_rate
    float(ar[3]) * (old_interval[1] - old_interval[0]) / (int(ar[2]) - int(ar[1]))
ZeroDivisionError: float division by zero
```


Discussion / hints:
Possibly related to #595.
`liftOver` is outputting 0-length bed intervals.
```
awk '$2==$3' genetic_map_PonAbe3_1.txt.tmp2 |head
chr1    344952  344952  0.261933        int2723
chr1    963630  963630  0.004067        int6586
chr1    2098258 2098258 11.502701       int13362
chr1    4132776 4132776 1.317146        int28819
chr1    4136651 4136651 0.373152        int28849
chr1    4138772 4138772 0.337640        int28865
chr1    4157123 4157123 0.636030        int28996
chr1    4176743 4176743 3.570152        int29126
chr1    4192872 4192872 0.731573        int29296
chr1    8095452 8095452 0.019146        int60004
```
The `hapmap_to_bed()` function in `maintenance/liftOver.py` is producing 0-length intervals, which are then given to `liftOver`. I guess there's an off-by-one error here, but the code is hard to follow.
https://github.com/popsim-consortium/stdpopsim/blob/5c149f9914823f3afef3f4ada8cf4c320edb5bbd/maintenance/liftOver.py#L44-L97
It's a problem with the recombination map, which itself contains 0-length intervals. How is one to interpret this (3rd column is the rate (cM/Mb), the 4th column is the genetic position (cM))?
```
awk '$2==524260' ~/.cache/stdpopsim/genetic_maps/PonAbe/NaterPA_PonAbe2/Nater_et_al_PA_chr1_PonAbe2.txt
chr1    524260  0.261933        0.246834
chr1    524260  0.322591        0.246834
```
Um wow that's weird. I think we just have to take it at face value, though, and asusme that that recombination rate applies to zero base pairs? (It's probably a bug generating the map somehow, but we're not trying to debug that?)

Fix patch:
diff --git a/stdpopsim/catalog/PonAbe/demographic_models.py b/stdpopsim/catalog/PonAbe/demographic_models.py
index 675d06bca..32c19b9ce 100644
--- a/stdpopsim/catalog/PonAbe/demographic_models.py
+++ b/stdpopsim/catalog/PonAbe/demographic_models.py
@@ -39,7 +39,9 @@ def _orangutan():
     # time of split
     T_split_years = 403149
     # get split time in units of generations
-    generation_time = _species.generation_time
+    # generation time used by Locke et al. without further citation
+    generation_time = 20
+    # mutation rate used by Locke et al. without further citation
     mutation_rate = 2e-8
 
     T_split = T_split_years / generation_time
diff --git a/stdpopsim/catalog/PonAbe/genetic_maps.py b/stdpopsim/catalog/PonAbe/genetic_maps.py
index e077e7679..e10460c7c 100644
--- a/stdpopsim/catalog/PonAbe/genetic_maps.py
+++ b/stdpopsim/catalog/PonAbe/genetic_maps.py
@@ -13,42 +13,45 @@
 # Recombination maps from Nater et al. were converted from rho/kbp to cM using
 # Watterson's estimator of theta to estimate Ne = 41,000 (Sumatra) and
 # Ne = 27,000 (Borneo). See supporting information in Nater et al. for details.
+# Maps have been lifted over to PonAbe3 (aka Susie_PABv2).
 
 
 _gm_pa = stdpopsim.GeneticMap(
     species=_species,
-    id="NaterPA_PonAbe2",
+    id="NaterPA_PonAbe3",
     description="From Nater et al. (2017) for Pongo abelii",
     long_description="""
         This genetic map is from the Nater et al. (2017) study, inferred using
         LDhat from n=15 whole-genome sequenced Sumatran orangutan individuals.
         See https://doi.org/10.1016/j.cub.2017.09.047 for more details.
+        Lifted over from assembly PonAbe2 (as used in Nater et al.) to PonAbe3.
         """,
     url=(
         "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonAbe/"
-        "NaterPA_PonAbe2.tar.gz"
+        "NaterPA_PonAbe3.tar.gz"
     ),
-    sha256="33b0162a6d945dd341bd1086d213ad1bc16949c9210d3b49d92692fd7f831ace",
-    file_pattern="Nater_et_al_PA_chr{id}_PonAbe2.txt",
+    sha256="80470e8d0a237f9cc938cac4595b13a6511b6988706a86db905c43e4619dfe37",
+    file_pattern="genetic_map_PonAbe3_{id}.txt",
     citations=[_nater2017.because(stdpopsim.CiteReason.GEN_MAP)],
 )
 _species.add_genetic_map(_gm_pa)
 
 _gm_pp = stdpopsim.GeneticMap(
     species=_species,
-    id="NaterPP_PonAbe2",
+    id="NaterPP_PonAbe3",
     description="From Nater et al. (2017) for Pongo pygmaeus",
     long_description="""
         This genetic map is from the Nater et al. (2017) study, inferred using
         LDhat from n=20 whole-genome sequenced Bornean orangutan individuals.
         See https://doi.org/10.1016/j.cub.2017.09.047 for more details.
+        Lifted over from assembly PonAbe2 (as used in Nater et al.) to PonAbe3.
         """,
     url=(
-        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonPyg/"
-        "NaterPP_PonAbe2.tar.gz"
+        "https://stdpopsim.s3-us-west-2.amazonaws.com/genetic_maps/PonAbe/"
+        "NaterPP_PonAbe3.tar.gz"
     ),
-    sha256="f4858b7efe15abe28b9367e7e9dc16614dc614df9326f9eedbcd943ca075bed7",
-    file_pattern="Nater_et_al_PP_chr{id}_PonAbe2.txt",
+    sha256="ccdaf29c3a7e515e137bf16305da6da3f48f10d70998588335c650b01a622e4f",
+    file_pattern="genetic_map_PonAbe3_{id}.txt",
     citations=[_nater2017.because(stdpopsim.CiteReason.GEN_MAP)],
 )
 _species.add_genetic_map(_gm_pp)
diff --git a/stdpopsim/catalog/PonAbe/genome_data.py b/stdpopsim/catalog/PonAbe/genome_data.py
index dd2a39801..a329cef0a 100644
--- a/stdpopsim/catalog/PonAbe/genome_data.py
+++ b/stdpopsim/catalog/PonAbe/genome_data.py
@@ -1,32 +1,36 @@
-# File autogenerated from Ensembl REST API. Do not edit.
+# File autogenerated from Ensembl REST API.
+# WARNING: autogenerated file has been edited!
 data = {
-    "assembly_accession": None,
-    "assembly_name": "PPYG2",
+    "assembly_accession": "GCA_002880775.3",
+    "assembly_name": "ponAbe3",  # Also known as Susie_PABv2.
     "chromosomes": {
-        "1": {"length": 229942017, "synonyms": []},
-        "2a": {"length": 113028656, "synonyms": []},
-        "2b": {"length": 135000294, "synonyms": []},
-        "3": {"length": 202140232, "synonyms": []},
-        "4": {"length": 198332218, "synonyms": []},
-        "5": {"length": 183952662, "synonyms": []},
-        "6": {"length": 174210431, "synonyms": []},
-        "7": {"length": 157549271, "synonyms": []},
-        "8": {"length": 153482349, "synonyms": []},
-        "9": {"length": 135191526, "synonyms": []},
-        "10": {"length": 133410057, "synonyms": []},
-        "11": {"length": 132107971, "synonyms": []},
-        "12": {"length": 136387465, "synonyms": []},
-        "13": {"length": 117095149, "synonyms": []},
-        "14": {"length": 108868599, "synonyms": []},
-        "15": {"length": 99152023, "synonyms": []},
-        "16": {"length": 77800216, "synonyms": []},
-        "17": {"length": 73212453, "synonyms": []},
-        "18": {"length": 94050890, "synonyms": []},
-        "19": {"length": 60714840, "synonyms": []},
-        "20": {"length": 62736349, "synonyms": []},
-        "21": {"length": 48394510, "synonyms": []},
-        "22": {"length": 46535552, "synonyms": []},
-        "X": {"length": 156195299, "synonyms": []},
+        "1": {"length": 227913704, "synonyms": []},
+        # Capitalisation of 2A/2B changed from the previous assembly,
+        # so the old names are here as synonyms.
+        "2A": {"length": 109511694, "synonyms": ["2a", "chr2a"]},
+        "2B": {"length": 129937803, "synonyms": ["2b", "chr2b"]},
+        "3": {"length": 193656255, "synonyms": []},
+        "4": {"length": 189387572, "synonyms": []},
+        "5": {"length": 179185813, "synonyms": []},
+        "6": {"length": 169501136, "synonyms": []},
+        "7": {"length": 145408105, "synonyms": []},
+        "8": {"length": 144036388, "synonyms": []},
+        "9": {"length": 112206110, "synonyms": []},
+        "10": {"length": 132178492, "synonyms": []},
+        "11": {"length": 128122151, "synonyms": []},
+        "12": {"length": 132184051, "synonyms": []},
+        "13": {"length": 98475126, "synonyms": []},
+        "14": {"length": 88963417, "synonyms": []},
+        "15": {"length": 82547911, "synonyms": []},
+        "16": {"length": 68237989, "synonyms": []},
+        "17": {"length": 75914007, "synonyms": []},
+        "18": {"length": 75923960, "synonyms": []},
+        "19": {"length": 57575784, "synonyms": []},
+        "20": {"length": 60841859, "synonyms": []},
+        "21": {"length": 34683425, "synonyms": []},
+        "22": {"length": 35308119, "synonyms": []},
+        "X": {"length": 151242693, "synonyms": []},
+        # Mitochondria absent in ponAbe3, so length taken from ponAbe2.
         "MT": {"length": 16499, "synonyms": []},
     },
 }
diff --git a/stdpopsim/catalog/PonAbe/species.py b/stdpopsim/catalog/PonAbe/species.py
index 16201c919..7987b9482 100644
--- a/stdpopsim/catalog/PonAbe/species.py
+++ b/stdpopsim/catalog/PonAbe/species.py
@@ -3,43 +3,55 @@
 
 
 # Average recombination rates per chromosome taken from the Pongo abelii
-# recombination map inferred in Nater et al (doi:10.1016/j.cub.2017.09.047)
-# For chromosome X, we use the rate assumed in Locke et al.
+# recombination map inferred in Nater et al (doi:10.1016/j.cub.2017.09.047),
+# then lifted over from PonAbe2 to PonAbe3. PonAbe3 is nearly 400Mb shorter
+# so this makes a noticable difference to the calculated averages.
+# For chromosome X, we use the weighted genome-wide average.
 _recombination_rate_data = {
-    "1": 5.19e-9,
-    "2a": 5.43e-9,
-    "2b": 5.38e-9,
-    "3": 5.36e-9,
-    "4": 5.43e-9,
-    "5": 5.20e-9,
-    "6": 5.22e-9,
-    "7": 5.73e-9,
-    "8": 5.67e-9,
-    "9": 5.34e-9,
-    "10": 5.91e-9,
-    "11": 5.29e-9,
-    "12": 5.44e-9,
-    "13": 4.91e-9,
-    "14": 4.70e-9,
-    "15": 4.82e-9,
-    "16": 6.12e-9,
-    "17": 7.26e-9,
-    "18": 4.57e-9,
-    "19": 7.56e-9,
-    "20": 5.83e-9,
-    "21": 4.98e-9,
-    "22": 6.03e-9,
-    "X": 9.50e-9,
+    "1": 5.7209738520466636e-09,
+    "2A": 6.746942995474612e-09,
+    "2B": 6.088330307532069e-09,
+    "3": 5.776152240522997e-09,
+    "4": 6.061756714136327e-09,
+    "5": 5.8783946576216106e-09,
+    "6": 5.896045046261895e-09,
+    "7": 6.82979118342936e-09,
+    "8": 6.343129574251662e-09,
+    "9": 7.047556965001838e-09,
+    "10": 6.948119058906713e-09,
+    "11": 5.791845481175716e-09,
+    "12": 6.043159266363315e-09,
+    "13": 6.4081697383466255e-09,
+    "14": 6.110753350074597e-09,
+    "15": 6.631500663980084e-09,
+    "16": 6.7606578552415915e-09,
+    "17": 8.01982524161532e-09,
+    "18": 6.400648881415157e-09,
+    "19": 9.096468143371218e-09,
+    "20": 6.314242320359848e-09,
+    "21": 7.588672333415156e-09,
+    "22": 9.836284925758108e-09,
+    "X": 6.402884398829772e-09,
     "MT": 0,
 }
 
+# High-resolution comparative analysis of great ape genomes
+_kronenberg2018 = stdpopsim.Citation(
+    author="Kronenberg et al.",
+    year=2018,
+    doi="https://doi.org/10.1126/science.aar6343",
+    reasons={stdpopsim.CiteReason.ASSEMBLY},
+)
+
+# Comparative and demographic analysis of orang-utan genomes
 _locke2011 = stdpopsim.Citation(
     author="Locke et al.",
     year=2011,
     doi="http://doi.org/10.1038/nature09687",
-    reasons={stdpopsim.CiteReason.GEN_TIME, stdpopsim.CiteReason.POP_SIZE},
+    reasons={stdpopsim.CiteReason.POP_SIZE},
 )
 
+# Morphometric, Behavioral, and Genomic Evidence for a New Orangutan Species
 _nater2017 = stdpopsim.Citation(
     author="Nater et al.",
     year=2017,
@@ -47,6 +59,16 @@
     reasons={stdpopsim.CiteReason.MUT_RATE, stdpopsim.CiteReason.REC_RATE},
 )
 
+# Orangutans: Geographic Variation in Behavioral Ecology and Conservation
+# CHAPTER 5 Orangutan life history variation
+_wich2008 = stdpopsim.Citation(
+    author="Wich et al.",
+    year=2008,
+    doi="https://doi.org/10.1093/acprof:oso/9780199213276.003.0005",
+    reasons={stdpopsim.CiteReason.GEN_TIME},
+)
+
+
 _chromosomes = []
 for name, data in genome_data.data["chromosomes"].items():
     _chromosomes.append(
@@ -65,7 +87,7 @@
     chromosomes=_chromosomes,
     assembly_name=genome_data.data["assembly_name"],
     assembly_accession=genome_data.data["assembly_accession"],
-    citations=[_nater2017],
+    citations=[_kronenberg2018],
 )
 stdpopsim.utils.append_common_synonyms(_genome)
 
@@ -75,11 +97,11 @@
     name="Pongo abelii",
     common_name="Sumatran orangutan",
     genome=_genome,
-    # generation time used by Locke et al. without further citation
-    generation_time=20,
+    # generation time used by Nater et al., citing Wich et al.
+    generation_time=25,
     # Locke et al. inferred ancestral Ne
     population_size=1.79e4,
-    citations=[_locke2011],
+    citations=[_locke2011, _wich2008],
 )
 
 stdpopsim.register_species(_species)


Test patch:
diff --git a/tests/test_PonAbe.py b/tests/test_PonAbe.py
index 162b6b3aa..3faa4d6b1 100644
--- a/tests/test_PonAbe.py
+++ b/tests/test_PonAbe.py
@@ -1,22 +1,113 @@
+import pytest
+
 import stdpopsim
 from tests import test_species
 
 
-class TestSpecies(test_species.SpeciesTestBase):
+class TestSpeciesData(test_species.SpeciesTestBase):
+
     species = stdpopsim.get_species("PonAbe")
 
-    def test_basic_attributes(self):
+    def test_ensembl_id(self):
+        assert self.species.ensembl_id == "pongo_abelii"
+
+    def test_name(self):
+        assert self.species.name == "Pongo abelii"
+
+    def test_common_name(self):
+        assert self.species.common_name == "Sumatran orangutan"
+
+    def test_qc_population_size(self):
         assert self.species.population_size == 1.79 * 10**4
-        assert self.species.generation_time == 20
+
+    def test_qc_generation_time(self):
+        assert self.species.generation_time == 25
 
 
-class TestGenome(test_species.GenomeTestBase):
-    """
-    Tests for the Pongo abelii genome.
-    """
+class TestGenomeData(test_species.GenomeTestBase):
 
     genome = stdpopsim.get_species("PonAbe").genome
 
-    def test_basic_attributes(self):
-        # 23 + X + MT
-        assert len(self.genome.chromosomes) == 25
+    # Average rates obtained from the NaterPA_PonAbe3 maps.
+    # Loading all maps takes ~60 seconds, so we just hardcode the values.
+    # This is how they were obtained:
+    # # maps = {
+    # #     chrom.id: gmap.get_chromosome_map(chrom.id)
+    # #     for chrom in species.genome.chromosomes
+    # #     if chrom.id not in ("X", "MT")
+    # # }
+    # # mean_rates = {k: v.mean_rate for k, v in maps.items()}
+    @pytest.mark.parametrize(
+        ["name", "rate"],
+        {
+            "1": 5.7209738520466636e-09,
+            "2A": 6.746942995474612e-09,
+            "2B": 6.088330307532069e-09,
+            "3": 5.776152240522997e-09,
+            "4": 6.061756714136327e-09,
+            "5": 5.8783946576216106e-09,
+            "6": 5.896045046261895e-09,
+            "7": 6.82979118342936e-09,
+            "8": 6.343129574251662e-09,
+            "9": 7.047556965001838e-09,
+            "10": 6.948119058906713e-09,
+            "11": 5.791845481175716e-09,
+            "12": 6.043159266363315e-09,
+            "13": 6.4081697383466255e-09,
+            "14": 6.110753350074597e-09,
+            "15": 6.631500663980084e-09,
+            "16": 6.7606578552415915e-09,
+            "17": 8.01982524161532e-09,
+            "18": 6.400648881415157e-09,
+            "19": 9.096468143371218e-09,
+            "20": 6.314242320359848e-09,
+            "21": 7.588672333415156e-09,
+            "22": 9.836284925758108e-09,
+            # X isn't in the map, so this is the weighted genome-wide average.
+            # # np.average(
+            # #    [m.mean_rate for m in maps.values()],
+            # #    weights=[m.sequence_length for m in maps.values()],
+            # # )
+            "X": 6.402884398829772e-09,
+            "MT": 0,
+        }.items(),
+    )
+    def test_recombination_rate(self, name, rate):
+        assert rate == pytest.approx(
+            self.genome.get_chromosome(name).recombination_rate
+        )
+
+    _genome_mutation_rate = 1.5e-8
+
+    @pytest.mark.parametrize(
+        ["name", "rate"],
+        {
+            "1": _genome_mutation_rate,
+            "2a": _genome_mutation_rate,
+            "2b": _genome_mutation_rate,
+            "3": _genome_mutation_rate,
+            "4": _genome_mutation_rate,
+            "5": _genome_mutation_rate,
+            "6": _genome_mutation_rate,
+            "7": _genome_mutation_rate,
+            "8": _genome_mutation_rate,
+            "9": _genome_mutation_rate,
+            "10": _genome_mutation_rate,
+            "11": _genome_mutation_rate,
+            "12": _genome_mutation_rate,
+            "13": _genome_mutation_rate,
+            "14": _genome_mutation_rate,
+            "15": _genome_mutation_rate,
+            "16": _genome_mutation_rate,
+            "17": _genome_mutation_rate,
+            "18": _genome_mutation_rate,
+            "19": _genome_mutation_rate,
+            "20": _genome_mutation_rate,
+            "21": _genome_mutation_rate,
+            "22": _genome_mutation_rate,
+            "X": _genome_mutation_rate,
+            "MT": _genome_mutation_rate,  # XXX: this is almost certainly wrong!
+        }.items(),
+    )
+    def test_mutation_rate(self, name, rate):
+        assert rate == pytest.approx(self.genome.get_chromosome(name).mutation_rate)
diff --git a/tests/test_genetic_maps.py b/tests/test_genetic_maps.py
index 0b3c7331f..e1cf7e526 100644
--- a/tests/test_genetic_maps.py
+++ b/tests/test_genetic_maps.py
@@ -245,10 +245,6 @@ def test_unknown_chromosome(self):
     def test_one_chrom_from_each_map(self):
         for gm in stdpopsim.all_genetic_maps():
             species = gm.species
-            if gm.id in ("NaterPA_PonAbe2", "NaterPP_PonAbe2"):
-                # XXX: these maps are currently broken:
-                #   ValueError: The last entry in the 'rate' column must be zero
-                continue
             # Just load the first chromosome in the list.
             # There's no requirement that any given chromosome is actually
             # in the map, and we don't have a direct way to check