Skip to content

Commit eb8cdb3

Browse files
committed
Flesh out the basic workbooks with coalescence rate calcs and trees
1 parent b042c0e commit eb8cdb3

File tree

4 files changed

+503
-63
lines changed

4 files changed

+503
-63
lines changed

content/data/demo.arg

-373 KB
Binary file not shown.

content/data/demo.py

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
# demo.tree file created using the following script
2+
3+
import tskit
4+
import msprime
5+
import collections
6+
7+
from matplotlib import pyplot as plt
8+
9+
seq_length = 1e6
10+
recent_generations = 200
11+
rho = 1e-7
12+
tree = "((A:250,B:250):150,(C:210,D:210):190);"
13+
initial_size = collections.defaultdict(lambda: 500)
14+
initial_size.update({"A": 500, "B": 500, "C": 300, "D": 1000})
15+
samples = {"A": 10, "B": 10, "C": 10, "D": 10}
16+
17+
ts_store = []
18+
for name, Ne in initial_size.items():
19+
demography = msprime.Demography()
20+
demography.add_population(name=name, initial_size=Ne)
21+
# Place a selective sweep in population A
22+
if name == "A":
23+
p = 1 / (2 * Ne)
24+
sweep_model = msprime.SweepGenicSelection(
25+
position=seq_length//3,
26+
s=0.1,
27+
start_frequency=p,
28+
end_frequency=1 - p,
29+
dt=1 / (40 * Ne),
30+
)
31+
ts_store.append(msprime.sim_ancestry(
32+
samples[name],
33+
model=(sweep_model, msprime.StandardCoalescent()) if name == "A" else msprime.StandardCoalescent(),
34+
demography=demography,
35+
recombination_rate=rho,
36+
sequence_length=seq_length,
37+
end_time=recent_generations,
38+
random_seed=123,
39+
))
40+
tables = ts_store[0].dump_tables()
41+
for ts in ts_store[1:]:
42+
tables.union(ts.dump_tables(), node_mapping=[tskit.NULL] * ts.num_nodes)
43+
tables.provenances.clear()
44+
45+
demography = msprime.Demography.from_species_tree(tree, initial_size)
46+
ts = msprime.sim_ancestry(initial_state=tables.tree_sequence(), recombination_rate=rho, sequence_length=seq_length, demography=demography, random_seed=123).simplify()
47+
ts = msprime.sim_mutations(ts, rate=5e-8, random_seed=123)
48+
ts.dump("data/demo.trees")

content/data/demo.trees

100 KB
Binary file not shown.

0 commit comments

Comments
 (0)