diff --git a/stdpopsim/qc/MusMus.py b/stdpopsim/qc/MusMus.py new file mode 100644 index 000000000..10bbe85f0 --- /dev/null +++ b/stdpopsim/qc/MusMus.py @@ -0,0 +1,288 @@ +import msprime +import numpy as np +import stdpopsim + +_species = stdpopsim.get_species("MusMus") + +# gen time and mut rate for all three models +spec_generation_time = 1 +spec_mutation_rate = 5.7e-9 + +# Parameters are from Fujiwara et al. 2022 (Figure 3) +# The values themselves were provided by the authors +# directly to Peter Fields who implemented the model. +# The values were visually compared to the curves, +# and after a few tweaks, they were confirmed. +# Here, we use the same values in a separate +# implementation of the three demographic models + + +def QC_DomesticusEurope(): + # Domesticus model (blue curve in Fig 3 of + # Fujiwara et al. (2022) + id = "QC-DomesticusEurope_1F22" + pop_id = "M_musculus_domesticus" + time_n_size = np.array( + [ + (0, 2040), + (83, 3844), + (180, 90428), + (291, 145603), + (420, 111242), + (570, 115399), + (743, 147212), + (943, 159142), + (1175, 136620), + (1443, 97250), + (1754, 58488), + (2114, 33028), + (2530, 18939), + (3012, 11758), + (3570, 8463), + (4216, 7480), + (4964, 8332), + (5829, 11240), + (6831, 16490), + (7991, 23419), + (9334, 29931), + (10889, 34163), + (12688, 36886), + (14772, 41195), + (17183, 50557), + (19975, 67337), + (23207, 90926), + (26948, 115426), + (31279, 131016), + (36292, 132063), + (42096, 121751), + (48815, 107067), + (56592, 93046), + (65596, 81892), + (76019, 74185), + (88084, 69939), + (102052, 69317), + (118221, 73097), + (136938, 82953), + (158606, 101471), + (183689, 131392), + (212726, 173264), + (246340, 222951), + (285254, 271935), + (330300, 309961), + (382446, 327217), + (442812, 316861), + (512693, 279833), + (593589, 227037), + (687237, 173594), + (795646, 131050), + (921140, 98811), + (1066418, 98811), + (1234595, 133912), + (1429281, 133912), + (1654653, 133912), + ] + ) + + model = msprime.Demography() + model.add_population( + name=pop_id, description=pop_id, initial_size=time_n_size[0][1] + ) + for j in range(1, time_n_size.shape[0]): + time, size = time_n_size[j, :] + model.add_population_parameters_change( + time, + initial_size=size, + population=pop_id, + ) + + return stdpopsim.DemographicModel( + id=id, + description=id, + long_description=id, + generation_time=spec_generation_time, + mutation_rate=spec_mutation_rate, + model=model, + ) + + +def QC_MusculusKorea(): + # Musculus Korean population (red curve in Fig 3 of Fujiwara et al. (2022) + id = "QC-MusculusKorea_1F22" + pop_id = "M_musculus_musculus" + time_n_size = np.array( + [ + (0, 179912), + (35, 8931), + (76, 8035), + (123, 9029), + (177, 9960), + (240, 12104), + (313, 16254), + (398, 25527), + (495, 42715), + (609, 61935), + (740, 68111), + (891, 55959), + (1067, 36220), + (1270, 20382), + (1505, 11222), + (1778, 6695), + (2093, 4605), + (2458, 3751), + (2881, 3643), + (3370, 4177), + (3936, 5506), + (4591, 7990), + (5350, 12072), + (6229, 17741), + (7246, 23546), + (8423, 26648), + (9785, 25399), + (11363, 21219), + (13189, 16747), + (15303, 13588), + (17750, 12259), + (20583, 13023), + (23863, 16339), + (27659, 22556), + (32054, 30806), + (37142, 38441), + (43031, 42857), + (49849, 43874), + (57741, 43467), + (66878, 43933), + (77455, 47001), + (89698, 54304), + (103872, 67725), + (120280, 88494), + (139274, 116547), + (161262, 151909), + (186716, 194969), + (216182, 245823), + (250293, 302950), + (289781, 359368), + (335491, 400867), + (388409, 407105), + (449667, 407105), + (520579, 152757), + (602670, 152757), + (697702, 152757), + ] + ) + + model = msprime.Demography() + model.add_population( + name=pop_id, description=pop_id, initial_size=time_n_size[0][1] + ) + for j in range(1, time_n_size.shape[0]): + time, size = time_n_size[j, :] + model.add_population_parameters_change( + time, + initial_size=size, + population=pop_id, + ) + + return stdpopsim.DemographicModel( + id=id, + description=id, + long_description=id, + generation_time=spec_generation_time, + mutation_rate=spec_mutation_rate, + model=model, + ) + + +def QC_CastaneusIndia(): + # Castaneus Indian model (green curve in Fig 3 of Fujiwara et al. (2022) + id = "QC-CastaneusIndia_1F22" + pop_id = "M_musculus_castaneus" + time_n_size = np.array( + [ + (0, 64853), + (887, 5064712), + (1886, 938111), + (3011, 291323), + (4279, 141377), + (5709, 86633), + (7319, 60911), + (9133, 47736), + (11178, 41845), + (13481, 41297), + (16077, 45372), + (19002, 53747), + (22297, 66260), + (26011, 83216), + (30195, 105533), + (34909, 134861), + (40221, 173419), + (46207, 222464), + (52951, 279816), + (60551, 338426), + (69114, 388298), + (78762, 420389), + (89634, 429975), + (101884, 418882), + (115687, 394606), + (131239, 365787), + (148764, 338818), + (168510, 317235), + (190760, 302369), + (215830, 294117), + (244079, 291560), + (275907, 293233), + (311772, 297138), + (352184, 300840), + (397719, 301783), + (449026, 297772), + (506839, 287572), + (571979, 271326), + (645379, 250544), + (728084, 227605), + (821274, 205098), + (926277, 185260), + (1044593, 169657), + (1177907, 159139), + (1328123, 154098), + (1497384, 155027), + (1688102, 163295), + (1903000, 182054), + (2145140, 217017), + (2417965, 275733), + (2725404, 360433), + (3071789, 463464), + (3462105, 463464), + (3901912, 344802), + (4397456, 344802), + (4955842, 344802), + ] + ) + + model = msprime.Demography() + model.add_population( + name=pop_id, description=pop_id, initial_size=time_n_size[0][1] + ) + for j in range(1, time_n_size.shape[0]): + time, size = time_n_size[j, :] + model.add_population_parameters_change( + time, + initial_size=size, + population=pop_id, + ) + + return stdpopsim.DemographicModel( + id=id, + description=id, + long_description=id, + generation_time=spec_generation_time, + mutation_rate=spec_mutation_rate, + model=model, + ) + + +_species.get_demographic_model("DomesticusEurope_1F22").register_qc( + QC_DomesticusEurope() +) + +_species.get_demographic_model("MusculusKorea_1F22").register_qc(QC_MusculusKorea()) + +_species.get_demographic_model("CastaneusIndia_1F22").register_qc(QC_CastaneusIndia()) diff --git a/stdpopsim/qc/__init__.py b/stdpopsim/qc/__init__.py index 847f9745b..20b4013b2 100644 --- a/stdpopsim/qc/__init__.py +++ b/stdpopsim/qc/__init__.py @@ -8,3 +8,4 @@ from . import PapAnu # NOQA from . import PanTro # NOQA from . import OrySat # NOQA +from . import MusMus # NOQA