forked from popsim-consortium/stdpopsim
-
Notifications
You must be signed in to change notification settings - Fork 0
/
dmel_example.py
67 lines (51 loc) · 1.84 KB
/
dmel_example.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
"""
Example of using the stdpopsim library with msprime.
"""
import msprime
import stdpopsim
from stdpopsim import drosophila_melanogaster
chrom = drosophila_melanogaster.genome.chromosomes["chr2R"]
recomb_map = chrom.recombination_map()
print("Testing Li and Stephan (2006) model")
model = drosophila_melanogaster.LiStephanTwoPopulation()
model.debug()
samples = [msprime.Sample(population=0, time=0),
msprime.Sample(population=0, time=0),
msprime.Sample(population=1, time=0),
msprime.Sample(population=1, time=0)]
ts = msprime.simulate(
samples=samples,
# recombination rate placeholder for quick runtime
recombination_rate=1e-09,
mutation_rate=chrom.default_mutation_rate,
**model.asdict())
print("simulated:", ts.num_trees, ts.num_sites)
print("====================\n")
print("Testing Sheehan and Song (2016) model")
model = drosophila_melanogaster.SheehanSongThreeEpoch()
model.debug()
samples = [msprime.Sample(population=0, time=0),
msprime.Sample(population=0, time=0)
]
ts = msprime.simulate(
samples=samples,
# recombination rate placeholder for quick runtime
recombination_rate=1e-09,
mutation_rate=chrom.default_mutation_rate,
**model.asdict())
print("simulated:", ts.num_trees, ts.num_sites)
print("====================\n")
print("Testing Sheehan and Song (2016) model on chrom2R map.")
print("This will take a while- perhaps 2 hours\n")
model = drosophila_melanogaster.SheehanSongThreeEpoch()
model.debug()
samples = [msprime.Sample(population=0, time=0),
msprime.Sample(population=0, time=0)
]
ts = msprime.simulate(
samples=samples,
# actual recombination map; slow
recombination_map=chrom.recombination_map(),
mutation_rate=chrom.default_mutation_rate,
**model.asdict())
print("simulated:", ts.num_trees, ts.num_sites)