Merge pull request #720 from grahamgower/chrom_length

Check the genetic map length matches the chromsome length.
popsim-consortium · Dec 17, 2020 · ec32ab9 · ec32ab9
2 parents 472f3f7 + 0a8c708
commit ec32ab9
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 14 deletions.
diff --git a/stdpopsim/genetic_maps.py b/stdpopsim/genetic_maps.py
@@ -100,14 +100,25 @@ def get_chromosome_map(self, id):
         # needs to be redownloaded.
         map_file = self.map_cache_dir / self.file_pattern.format(id=chrom.id)
         if map_file.exists():
-            ret = msprime.RecombinationMap.read_hapmap(str(map_file))
+            recomb_map = msprime.RecombinationMap.read_hapmap(str(map_file))
         else:
             warnings.warn(
-                "Warning: recombination map not found for chromosome: '{}'"
+                "Recombination map not found for chromosome: '{}'"
                 " on map: '{}', substituting a flat map with chromosome "
                 "recombination rate {}".format(id, self.id, chrom.recombination_rate)
             )
-            ret = msprime.RecombinationMap.uniform_map(
+            recomb_map = msprime.RecombinationMap.uniform_map(
                 chrom.length, chrom.recombination_rate
             )
-        return ret
+        map_length = recomb_map.get_sequence_length()
+        if map_length < chrom.length:
+            # Extend map to the end of the chromosome.
+            positions = recomb_map.get_positions() + [chrom.length]
+            rates = recomb_map.get_rates() + [0]
+            recomb_map = msprime.RecombinationMap(positions, rates)
+        elif map_length > chrom.length:
+            warnings.warn(
+                f"Recombination map has length {map_length}, which is longer than"
+                f" chomosome length {chrom.length}. The former will be used."
+            )
+        return recomb_map
diff --git a/tests/test_genetic_maps.py b/tests/test_genetic_maps.py
@@ -202,22 +202,35 @@ class TestGetChromosomeMap(tests.CacheReadingTest):
     Tests if we get chromosome maps using the HapMapII_GRCh37 human map.
     """
 
-    species = stdpopsim.get_species("HomSap")
-    genetic_map = species.get_genetic_map("HapMapII_GRCh37")
-
     def test_warning_from_no_mapped_chromosome(self):
-        chrom = self.species.genome.get_chromosome("chrY")
+        species = stdpopsim.get_species("HomSap")
+        genetic_map = species.get_genetic_map("HapMapII_GRCh37")
+        chrom = species.genome.get_chromosome("chrY")
         with self.assertWarns(Warning):
-            cm = self.genetic_map.get_chromosome_map(chrom.id)
-            self.assertIsInstance(cm, msprime.RecombinationMap)
-            self.assertEqual(chrom.length, cm.get_sequence_length())
+            cm = genetic_map.get_chromosome_map(chrom.id)
+        self.assertIsInstance(cm, msprime.RecombinationMap)
+        self.assertEqual(chrom.length, cm.get_sequence_length())
 
     def test_known_chromosome(self):
-        chrom = self.species.genome.get_chromosome("chr22")
-        cm = self.genetic_map.get_chromosome_map(chrom.id)
+        species = stdpopsim.get_species("CanFam")
+        genetic_map = species.get_genetic_map("Campbell2016_CanFam3_1")
+        chrom = species.genome.get_chromosome("1")
+        cm = genetic_map.get_chromosome_map(chrom.id)
+        self.assertIsInstance(cm, msprime.RecombinationMap)
+        self.assertEqual(chrom.length, cm.get_sequence_length())
+
+    def test_warning_for_long_recomb_map(self):
+        species = stdpopsim.get_species("HomSap")
+        genetic_map = species.get_genetic_map("HapMapII_GRCh37")
+        chrom = species.genome.get_chromosome("chr1")
+        with self.assertWarns(Warning):
+            cm = genetic_map.get_chromosome_map(chrom.id)
         self.assertIsInstance(cm, msprime.RecombinationMap)
+        self.assertLess(chrom.length, cm.get_sequence_length())
 
     def test_unknown_chromosome(self):
+        species = stdpopsim.get_species("HomSap")
+        genetic_map = species.get_genetic_map("HapMapII_GRCh37")
         for bad_chrom in ["", "ABD", None]:
             with self.assertRaises(ValueError):
-                self.genetic_map.get_chromosome_map(bad_chrom)
+                genetic_map.get_chromosome_map(bad_chrom)