add mus mus domesticus demographic model #1485
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,344 @@ | ||
| import msprime | ||
| import numpy as np | ||
|
|
||
| import stdpopsim | ||
|
|
||
| _species = stdpopsim.get_species("MusMus") | ||
|
|
||
|
|
||
| def _dom_1pop(): | ||
| # the size during the interval times[k] to times[k+1] = sizes[k] | ||
| times = np.array( | ||
| [ | ||
| 0, | ||
| 83, | ||
| 180, | ||
| 291, | ||
| 420, | ||
| 570, | ||
| 743, | ||
| 943, | ||
| 1175, | ||
| 1443, | ||
| 1754, | ||
| 2114, | ||
| 2530, | ||
| 3012, | ||
| 3570, | ||
| 4216, | ||
| 4964, | ||
| 5829, | ||
| 6831, | ||
| 7991, | ||
| 9334, | ||
| 10889, | ||
| 12688, | ||
| 14772, | ||
| 17183, | ||
| 19975, | ||
| 23207, | ||
| 26948, | ||
| 31279, | ||
| 36292, | ||
| 42096, | ||
| 48815, | ||
| 56592, | ||
| 65596, | ||
| 76019, | ||
| 88084, | ||
| 102052, | ||
| 118221, | ||
| 136938, | ||
| 158606, | ||
| 183689, | ||
| 212726, | ||
| 246340, | ||
| 285254, | ||
| 330300, | ||
| 382446, | ||
| 442812, | ||
| 512693, | ||
| 593589, | ||
| 687237, | ||
| 795646, | ||
| 921140, | ||
| 1066418, | ||
| 1234595, | ||
| 1429281, | ||
| 1654653, | ||
| 1915544, | ||
| ] | ||
| ) | ||
| sizes = np.array( | ||
| [ | ||
| 2040, | ||
| 2040, | ||
| 3844, | ||
| 90428, | ||
| 145603, | ||
| 111242, | ||
| 115399, | ||
| 147212, | ||
| 159142, | ||
| 136620, | ||
| 97250, | ||
| 58488, | ||
| 33028, | ||
| 18939, | ||
| 11758, | ||
| 8463, | ||
| 7480, | ||
| 8332, | ||
| 11240, | ||
| 16490, | ||
| 23419, | ||
| 29931, | ||
| 34163, | ||
| 36886, | ||
| 41195, | ||
| 50557, | ||
| 67337, | ||
| 90926, | ||
| 115426, | ||
| 131016, | ||
| 132063, | ||
| 121751, | ||
| 107067, | ||
| 93046, | ||
| 81892, | ||
| 74185, | ||
| 69939, | ||
| 69317, | ||
| 73097, | ||
| 82953, | ||
| 101471, | ||
| 131392, | ||
| 173264, | ||
| 222951, | ||
| 271935, | ||
| 309961, | ||
| 327217, | ||
| 316861, | ||
| 279833, | ||
| 227037, | ||
| 173594, | ||
| 131050, | ||
| 98811, | ||
| 98811, | ||
| 133912, | ||
| 133912, | ||
| 133912, | ||
| ] | ||
| ) | ||
|
|
||
| demographic_events = [] | ||
| for sz, t in zip(sizes, times): | ||
| demographic_events.append( | ||
| msprime.PopulationParametersChange(time=t, initial_size=sz, population_id=0) | ||
| ) | ||
|
|
||
| populations = [ | ||
| stdpopsim.Population( | ||
| id="M_musculus_domesticus", | ||
| description="Mus musculus domesticus German population", | ||
| ) | ||
| ] | ||
|
|
||
| return stdpopsim.DemographicModel( | ||
| id="M_musculus_domesticus_Europe", | ||
| description="M. musculus domesticus piecewise constant size", | ||
|
There was a problem hiding this comment. We have certain conventions for the demographic model ids. You can find them described here. According to these conventions, I would suggest something like DomesticusEurope_1F22 There was a problem hiding this comment. Done. I also added similar names for the other demographic models. |
||
| long_description=""" | ||
| This model comes from MSMC using four randomly sampled | ||
| individuals (DEU01,DEU03,DEU04,DEU06) from a German population. | ||
| The model is estimated with 57 time periods. | ||
| """, | ||
|
There was a problem hiding this comment. I would also explicitly mention Fig 3 of Fujiwara et al 2022 and say that the population sizes and time changes were supplied by the authors (maybe even mention a specific author?). There was a problem hiding this comment. I have added these in. |
||
| populations=populations, | ||
| citations=[ | ||
| stdpopsim.Citation( | ||
| author="Fujiwara et al.", | ||
| year=2022, | ||
| doi="https://doi.org/10.1093/gbe/evac068", | ||
| reasons={stdpopsim.CiteReason.DEM_MODEL}, | ||
| ) | ||
| ], | ||
| generation_time=1, | ||
| mutation_rate=5.7e-9, | ||
| demographic_events=demographic_events, | ||
| population_configurations=[ | ||
| msprime.PopulationConfiguration( | ||
| initial_size=sizes[0], metadata=populations[0].asdict() | ||
| ) | ||
| ], | ||
| ) | ||
|
|
||
|
|
||
| _species.add_demographic_model(_dom_1pop()) | ||
|
|
||
|
|
||
| def _mus_1pop(): | ||
| # the size during the interval times[k] to times[k+1] = sizes[k] | ||
| times = np.array( | ||
| [ | ||
| 0, | ||
| 35, | ||
| 76, | ||
| 123, | ||
| 177, | ||
| 240, | ||
| 313, | ||
| 398, | ||
| 495, | ||
| 609, | ||
| 740, | ||
| 891, | ||
| 1067, | ||
| 1270, | ||
| 1505, | ||
| 1778, | ||
| 2093, | ||
| 2458, | ||
| 2881, | ||
| 3370, | ||
| 3936, | ||
| 4591, | ||
| 5350, | ||
| 6229, | ||
| 7246, | ||
| 8423, | ||
| 9785, | ||
| 11363, | ||
| 13189, | ||
| 15303, | ||
| 17750, | ||
| 20583, | ||
| 23863, | ||
| 27659, | ||
| 32054, | ||
| 37142, | ||
| 43031, | ||
| 49849, | ||
| 57741, | ||
| 66878, | ||
| 77455, | ||
| 89698, | ||
| 103872, | ||
| 120280, | ||
| 139274, | ||
| 161262, | ||
| 186716, | ||
| 216182, | ||
| 250293, | ||
| 289781, | ||
| 335491, | ||
| 388409, | ||
| 449667, | ||
| 520579, | ||
| 602670, | ||
| 697702, | ||
| 807711, | ||
| ] | ||
| ) | ||
| sizes = np.array( | ||
| [ | ||
| 179912, | ||
| 179912, | ||
| 8931, | ||
| 8035, | ||
| 9029, | ||
| 9960, | ||
| 12104, | ||
| 16254, | ||
| 25527, | ||
| 42715, | ||
| 61935, | ||
| 68111, | ||
| 55959, | ||
| 36220, | ||
| 20382, | ||
| 11222, | ||
| 6695, | ||
| 4605, | ||
| 3751, | ||
| 3643, | ||
| 4177, | ||
| 5506, | ||
| 7990, | ||
| 12072, | ||
| 17741, | ||
| 23546, | ||
| 26648, | ||
| 25399, | ||
| 21219, | ||
| 16747, | ||
| 13588, | ||
| 12259, | ||
| 13023, | ||
| 16339, | ||
| 22556, | ||
| 30806, | ||
| 38441, | ||
| 42857, | ||
| 43874, | ||
| 43467, | ||
| 43933, | ||
| 47001, | ||
| 54304, | ||
| 67725, | ||
| 88494, | ||
| 116547, | ||
| 151909, | ||
| 194969, | ||
| 245823, | ||
| 302950, | ||
| 359368, | ||
| 400867, | ||
| 407105, | ||
| 407105, | ||
| 152757, | ||
| 152757, | ||
| 152757, | ||
| ] | ||
| ) | ||
|
|
||
| demographic_events = [] | ||
| for sz, t in zip(sizes, times): | ||
| demographic_events.append( | ||
| msprime.PopulationParametersChange(time=t, initial_size=sz, population_id=0) | ||
| ) | ||
|
|
||
| populations = [ | ||
| stdpopsim.Population( | ||
| id="M_musculus_musculus", | ||
| description="Mus musculus musculus Korean population", | ||
| ) | ||
| ] | ||
|
|
||
| return stdpopsim.DemographicModel( | ||
| id="M_musculus_musculus_East_Asia", | ||
| description="M. musculus musculus piecewise constant size", | ||
| long_description=""" | ||
| This model comes from MSMC using four randomly sampled | ||
| individuals (KOR01,KOR02,KOR03,KOR05) from a Korean population. | ||
| The model is estimated with 57 time periods. | ||
| """, | ||
|
There was a problem hiding this comment. See comments above about the model id and long description |
||
| populations=populations, | ||
| citations=[ | ||
| stdpopsim.Citation( | ||
| author="Fujiwara et al.", | ||
| year=2022, | ||
| doi="https://doi.org/10.1093/gbe/evac068", | ||
| reasons={stdpopsim.CiteReason.DEM_MODEL}, | ||
| ) | ||
| ], | ||
| generation_time=1, | ||
| mutation_rate=5.7e-9, | ||
| demographic_events=demographic_events, | ||
| population_configurations=[ | ||
| msprime.PopulationConfiguration( | ||
| initial_size=sizes[0], metadata=populations[0].asdict() | ||
| ) | ||
| ], | ||
| ) | ||
|
|
||
|
|
||
| _species.add_demographic_model(_mus_1pop()) | ||
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
check that there is no offset here. From examination of fig 3, it looks like between 300-400 generations ago the Ne was ~150K. The way the arrays are set up here, I think that Ne=145,603 is associated with the time range 420-570.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
When you say offset here do you mean something being added at the stage of plotting?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There appears to be an inconsistency between the plot in fig 3 of Fujiwara et al 2022 and your implementation. I don't know if the problem is in fig 3 or your implementation. If I understand the implementation correctly, then each Ne in your table is being associated with the wrong time range (shifted one range back in time). The 5th Ne (145603) is associated in your implementation with the 5th time interval (420 - 570 generations ago). However, in Fig. 3, it appears to be associated the time range somewhere between 300 to 400 generations ago, which fits your 4th time interval. If you remove the first element in the
sizesarray, this should fix it IMO.