Add ploidy as attribute of species; rework samples specification (#1361)

* Ploidy as species/chromosome attribute; use dict for samples in python API; use population:integer pairs for CLI

* Use custom deprecation warning for old samples API

docs/examples.py

@@ -13,7 +13,7 @@ def generic_models_example():
     species = stdpopsim.get_species("HomSap")
     contig = species.get_contig("chr22", length_multiplier=0.1)
     model = stdpopsim.PiecewiseConstantSize(species.population_size)
-    samples = model.get_samples(10)
+    samples = {"pop_0": 5}
     engine = stdpopsim.get_default_engine()
     ts = engine.simulate(model, contig, samples)
     print("num_trees =", ts.num_trees)

docs/selection_example.py

@@ -21,9 +21,7 @@ def adaptive_introgression(seed):
     species = stdpopsim.get_species("HomSap")
     model = species.get_demographic_model("PapuansOutOfAfrica_10J19")
     contig = species.get_contig("chr1", length_multiplier=0.001)
-    samples = model.get_samples(
-        100, 0, 0, 100, 2, 2  # YRI, CEU, CHB, Papuan, DenA, NeaA
-    )
+    samples = {"YRI": 50, "Papuan": 50, "DenA": 1, "NeaA": 1}
 
     # We need some demographic model parameters to set bounds on the timing
     # of random variables and extended_events (below).

stdpopsim/catalog/AedAeg/species.py

@@ -5,6 +5,15 @@
 # These are in Table 1 of Juneja et al:
 _recombination_rate = {"1": 0.306e-8, "2": 0.249e-8, "3": 0.291e-8, "MT": 0}
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "1": _species_ploidy,
+    "2": _species_ploidy,
+    "3": _species_ploidy,
+    "MT": 1,
+}
+
 _JunejaEtAl = stdpopsim.Citation(
     doi="https://doi.org/10.1371/journal.pntd.0002652",
     year=2014,
@@ -56,6 +65,7 @@
     genome_data.data,
     recombination_rate=_recombination_rate,
     mutation_rate=_mutation_rate,
+    ploidy=_ploidy,
     citations=[
         _NeneEtAl,
         _JunejaEtAl,
@@ -72,6 +82,7 @@
     common_name="Yellow fever mosquito",
     genome=_genome,
     generation_time=1 / 15,
+    ploidy=_species_ploidy,
     # the estimated population size today the modern Senegal forest population
     population_size=1e6,
     citations=[_CrawfordEtAl],

stdpopsim/catalog/AnaPla/species.py

@@ -148,6 +148,51 @@
     "40": _default_recombination_rate,
 }
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "1": _species_ploidy,
+    "2": _species_ploidy,
+    "3": _species_ploidy,
+    "4": _species_ploidy,
+    "5": _species_ploidy,
+    "6": _species_ploidy,
+    "7": _species_ploidy,
+    "8": _species_ploidy,
+    "9": _species_ploidy,
+    "10": _species_ploidy,
+    "11": _species_ploidy,
+    "12": _species_ploidy,
+    "13": _species_ploidy,
+    "14": _species_ploidy,
+    "15": _species_ploidy,
+    "16": _species_ploidy,
+    "17": _species_ploidy,
+    "18": _species_ploidy,
+    "19": _species_ploidy,
+    "20": _species_ploidy,
+    "21": _species_ploidy,
+    "22": _species_ploidy,
+    "23": _species_ploidy,
+    "24": _species_ploidy,
+    "25": _species_ploidy,
+    "26": _species_ploidy,
+    "27": _species_ploidy,
+    "28": _species_ploidy,
+    "29": _species_ploidy,
+    "30": _species_ploidy,
+    "Z": _species_ploidy,
+    "31": _species_ploidy,
+    "32": _species_ploidy,
+    "33": _species_ploidy,
+    "34": _species_ploidy,
+    "35": _species_ploidy,
+    "36": _species_ploidy,
+    "37": _species_ploidy,
+    "38": _species_ploidy,
+    "39": _species_ploidy,
+    "40": _species_ploidy,
+}
 
 # value used in Lavretsky et al 2020, as obtained for nuclear genes in other
 # ducks by Peters, Zhuravlev, Fefelov, Humphries, & Omland, 2008
@@ -203,6 +248,7 @@
     genome_data.data,
     recombination_rate=_recombination_rate,
     mutation_rate=_mutation_rate,
+    ploidy=_ploidy,
     citations=[
         _LavretskyEtAl2019,
         _HuangEtAl2006,
@@ -237,6 +283,7 @@
     # choosing Ne based on theta = 4 Ne u from Guo et al 2021
     # theta = 0.003 (Figure 1), u as above (the paper uses a rate from chicken)
     population_size=156000,
+    ploidy=_species_ploidy,
     citations=[
         _LavretskyEtAl2020,
         _GuoEtAl2020,

stdpopsim/catalog/AnoCar/species.py

@@ -42,6 +42,25 @@
     "MT": 0,
 }
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "1": _species_ploidy,
+    "2": _species_ploidy,
+    "3": _species_ploidy,
+    "4": _species_ploidy,
+    "5": _species_ploidy,
+    "6": _species_ploidy,
+    "LGa": _species_ploidy,
+    "LGb": _species_ploidy,
+    "LGc": _species_ploidy,
+    "LGd": _species_ploidy,
+    "LGf": _species_ploidy,
+    "LGg": _species_ploidy,
+    "LGh": _species_ploidy,
+    "MT": 1,
+}
+
 # Mutation rate
 _overall_rate = 2.1e-10
 
@@ -66,6 +85,7 @@
     genome_data.data,
     recombination_rate=_recombination_rate,
     mutation_rate=_mutation_rate,
+    ploidy=_ploidy,
     citations=[_BourgeoisEtAl],
 )
 stdpopsim.utils.append_common_synonyms(_genome)
@@ -84,6 +104,7 @@
     # poulation size caculated from theta caculations
     # theta = 4Neu, theta from table 1
     # Ne averaged across the 5 populations from BourgeoisEtAl
+    ploidy=_species_ploidy,
     citations=[_LovernEtAl, _BourgeoisEtAl],
 )
 

stdpopsim/catalog/AnoGam/species.py

@@ -56,6 +56,16 @@
     "Mt": 0,
 }
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "2L": _species_ploidy,
+    "2R": _species_ploidy,
+    "3L": _species_ploidy,
+    "3R": _species_ploidy,
+    "X": _species_ploidy,
+    "Mt": 1,
+}
 
 # the rate for DroMel as inferred by Keightley et al 2009,
 # which was used for the stairwayplot demographic inference
@@ -73,6 +83,7 @@
     genome_data.data,
     recombination_rate=_recombination_rate,
     mutation_rate=_mutation_rate,
+    ploidy=_ploidy,
     citations=[
         _SharakhovaEtAl,
         _Ag1000G,
@@ -92,6 +103,7 @@
     generation_time=1 / 11,
     # based on theta = 4 Ne u in Gabon population, rounded
     population_size=1e6,
+    ploidy=_species_ploidy,
     citations=[_Ag1000G],
 )
 

stdpopsim/catalog/ApiMel/species.py

@@ -30,6 +30,30 @@
     "CM009947.2": 0,
 }
 
+# TODO: This species is haplodiploid, but machinery is not in place to simulate
+# haplodiploidy.
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "CM009931.2": _species_ploidy,
+    "CM009932.2": _species_ploidy,
+    "CM009933.2": _species_ploidy,
+    "CM009934.2": _species_ploidy,
+    "CM009935.2": _species_ploidy,
+    "CM009936.2": _species_ploidy,
+    "CM009937.2": _species_ploidy,
+    "CM009938.2": _species_ploidy,
+    "CM009939.2": _species_ploidy,
+    "CM009940.2": _species_ploidy,
+    "CM009941.2": _species_ploidy,
+    "CM009942.2": _species_ploidy,
+    "CM009943.2": _species_ploidy,
+    "CM009944.2": _species_ploidy,
+    "CM009945.2": _species_ploidy,
+    "CM009946.2": _species_ploidy,
+    "CM009947.2": 1,  # Mt
+}
+
 
 _YangEtAl = stdpopsim.Citation(
     doi="https://doi.org/10.1038/nature14649",
@@ -79,6 +103,7 @@
     genome_data.data,
     recombination_rate=_recombination_rate,
     mutation_rate=_mutation_rate,
+    ploidy=_ploidy,
     citations=[
         _YangEtAl,
         _BeyeEtAl,
@@ -139,6 +164,7 @@ def add_if_unique(chrom, synonyms):
     genome=_genome,
     generation_time=2,
     population_size=2e05,
+    ploidy=_species_ploidy,
     citations=[
         _WallbergEtAl,
         _NelsonEtAl,

stdpopsim/catalog/AraTha/species.py

@@ -8,6 +8,18 @@
 _recombination_rate_data["Mt"] = 0
 _recombination_rate_data["Pt"] = 0  # JK Is this correct??
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "1": _species_ploidy,
+    "2": _species_ploidy,
+    "3": _species_ploidy,
+    "4": _species_ploidy,
+    "5": _species_ploidy,
+    "6": _species_ploidy,
+    "Mt": 1,
+    "Pt": 1,
+}
 
 # mutation rate from Ossowski 2010 Science
 # recombination value from Huber et al 2014 MBE
@@ -21,6 +33,7 @@
             synonyms=data["synonyms"],
             mutation_rate=7e-9,
             recombination_rate=_recombination_rate_data[name],
+            ploidy=_ploidy[name],
         )
     )
 
@@ -59,6 +72,7 @@
     genome=_genome,
     generation_time=1.0,
     population_size=10**4,
+    ploidy=_species_ploidy,
     citations=[
         stdpopsim.Citation(
             doi="https://doi.org/10.1890/0012-9658(2002)083[1006:GTINSO]2.0.CO;2",

stdpopsim/catalog/BosTau/species.py

@@ -64,6 +64,42 @@
 # Set some exceptions for non-recombining chrs.
 _recombination_rate_data["MT"] = 0
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "1": _species_ploidy,
+    "2": _species_ploidy,
+    "3": _species_ploidy,
+    "4": _species_ploidy,
+    "5": _species_ploidy,
+    "6": _species_ploidy,
+    "7": _species_ploidy,
+    "8": _species_ploidy,
+    "9": _species_ploidy,
+    "10": _species_ploidy,
+    "11": _species_ploidy,
+    "12": _species_ploidy,
+    "13": _species_ploidy,
+    "14": _species_ploidy,
+    "15": _species_ploidy,
+    "16": _species_ploidy,
+    "17": _species_ploidy,
+    "18": _species_ploidy,
+    "19": _species_ploidy,
+    "20": _species_ploidy,
+    "21": _species_ploidy,
+    "22": _species_ploidy,
+    "23": _species_ploidy,
+    "24": _species_ploidy,
+    "25": _species_ploidy,
+    "26": _species_ploidy,
+    "27": _species_ploidy,
+    "28": _species_ploidy,
+    "29": _species_ploidy,
+    "X": _species_ploidy,
+    "MT": 1,
+}
+
 _chromosomes = []
 for name, data in genome_data.data["chromosomes"].items():
     _chromosomes.append(
@@ -74,6 +110,7 @@
             # Harland et al. (2017), sex-averaged estimate per bp per generation.
             mutation_rate=1.2e-8,
             recombination_rate=_recombination_rate_data[name],
+            ploidy=_ploidy[name],
         )
     )
 
@@ -96,6 +133,7 @@
     genome=_genome,
     generation_time=5,
     population_size=62000,  # ancestral Ne in _MacLeodEtAl
+    ploidy=_species_ploidy,
     citations=[_MacLeodEtAl],
 )
 

stdpopsim/catalog/CaeEle/species.py

@@ -94,6 +94,18 @@
     "MtDNA": 1.05e-7,
 }
 
+# Generic and chromosome-specific ploidy
+_species_ploidy = 2
+_ploidy = {
+    "I": _species_ploidy,
+    "II": _species_ploidy,
+    "III": _species_ploidy,
+    "IV": _species_ploidy,
+    "V": _species_ploidy,
+    "X": _species_ploidy,
+    "MtDNA": 1,
+}
+
 _chromosomes = []
 for name, data in genome_data.data["chromosomes"].items():
     _chromosomes.append(
@@ -107,6 +119,7 @@
             # a mutation map.
             # mutation_rate=_mutation_rate_data[name],
             recombination_rate=_recombination_rate_data[name],
+            ploidy=_ploidy[name],
         )
     )
 
@@ -135,6 +148,7 @@
     population_size=10000,  # Lots of estimates available,
     # they are all over the place. Settling on 10,000.
     #    selfing_rate,=0.999,
+    ploidy=_species_ploidy,
     citations=[
         _FrezalAndFelix2015.because(stdpopsim.CiteReason.GEN_TIME),
         _BarriereAndFelix2005.because(stdpopsim.CiteReason.POP_SIZE),