swh:1:snp:300923221fcf626df34df8c763b7994a14d9c907
Tip revision: 772052e7ba53a55183f0bb9bd7ff9a4112e4e5c8 authored by Jerome Kelleher on 24 June 2016, 14:07:28 UTC
Added analytical verification for pi.
Added analytical verification for pi.
Tip revision: 772052e
examples.py
"""
The examples used in the tutorial section.
"""
from __future__ import print_function
from __future__ import division
import os
import sys
sys.path.insert(0, os.path.abspath('..'))
import math
import msprime
import numpy as np
import scipy.stats
import matplotlib
# Force matplotlib to not use any Xwindows backend.
matplotlib.use('Agg')
import matplotlib.pyplot as pyplot
import matplotlib.collections
def segregating_sites_example(n, theta, num_replicates):
S = np.zeros(num_replicates)
replicates = msprime.simulate(
sample_size=n,
mutation_rate=theta / 4,
num_replicates=num_replicates)
for j, tree_sequence in enumerate(replicates):
S[j] = tree_sequence.get_num_mutations()
S_mean_a = np.sum(1 / np.arange(1, n)) * theta
S_var_a = (
theta * np.sum(1 / np.arange(1, n)) +
theta**2 * np.sum(1 / np.arange(1, n)**2))
print(" mean variance")
print("Observed {}\t\t{}".format(np.mean(S), np.var(S)))
print("Analytical {:.5f}\t\t{:.5f}".format(S_mean_a, S_var_a))
def migration_example():
# M is the overall symmetric migration rate, and d is the number
# of demes.
M = 0.2
d = 3
# We rescale m into per-generation values for msprime.
m = M / (4 * (d - 1))
# Allocate the initial sample. Because we are interested in the
# between deme coalescence times, we choose one sample each
# from the first two demes.
population_configurations = [
msprime.PopulationConfiguration(sample_size=1),
msprime.PopulationConfiguration(sample_size=1),
msprime.PopulationConfiguration(sample_size=0)]
# Now we set up the migration matrix. Since this is a symmetric
# island model, we have the same rate of migration between all
# pairs of demes. Diagonal elements must be zero.
migration_matrix = [
[0, m, m],
[m, 0, m],
[m, m, 0]]
# We pass these values to the simulate function, and ask it
# to run the required number of replicates.
num_replicates = 1e6
replicates = msprime.simulate(
population_configurations=population_configurations,
migration_matrix=migration_matrix,
num_replicates=num_replicates)
# And then iterate over these replicates
T = np.zeros(num_replicates)
for i, tree_sequence in enumerate(replicates):
tree = next(tree_sequence.trees())
# Convert the TMRCA to coalecent units.
T[i] = tree.get_time(tree.get_root()) / 4
# Finally, calculate the analytical expectation and print
# out the results
analytical = d / 2 + (d - 1) / (2 * M)
print("Observed =", np.mean(T))
print("Predicted =", analytical)
def single_locus_example():
tree_sequence = msprime.simulate(sample_size=5, Ne=1000, random_seed=1)
tree = next(tree_sequence.trees())
print(tree)
tree.draw("_static/simple-tree.svg", show_times=True)
u = 0
while u != msprime.NULL_NODE:
print("node {}: time = {}".format(u, tree.get_time(u)))
u = tree.get_parent(u)
print(tree.get_branch_length(6))
print(tree.get_total_branch_length())
def multi_locus_example():
tree_sequence = msprime.simulate(
sample_size=5, Ne=1000, length=1e4, recombination_rate=2e-8,
random_seed=19)
j = 0
for tree in tree_sequence.trees():
print(tree.get_interval(), str(tree), sep="\t")
tree.draw("_static/simple-tree-sequence-{}.svg".format(j))
j += 1
def mutations_example():
tree_sequence = msprime.simulate(
sample_size=5, Ne=1000, length=1e4, recombination_rate=2e-8,
mutation_rate=2e-8, random_seed=19)
print("Total mutations = ", tree_sequence.get_num_mutations())
j = 0
for tree in tree_sequence.trees():
print(tree.get_interval(), list(tree.mutations()), sep="\t")
tree.draw("_static/mutations-tree-sequence-{}.svg".format(j))
j += 1
for tree in tree_sequence.trees():
for position, node in tree.mutations():
print("Mutation @ position {} has frequency {}".format(
position, tree.get_num_leaves(node) / tree.get_sample_size()))
def out_of_africa():
# First we set out the maximum likelihood values of the various parameters
# given in Table 1.
N_A = 7300
N_B = 2100
N_AF = 12300
N_EU0 = 1000
N_AS0 = 510
# Times are provided in years, so we convert into generations.
generation_time = 25
T_AF = 220e3 / generation_time
T_B = 140e3 / generation_time
T_EU_AS = 21.2e3 / generation_time
# We need to work out the starting population sizes based on the growth
# rates provided for these two populations
r_EU = 0.004
r_AS = 0.0055
N_EU = N_EU0 / math.exp(-r_EU * T_EU_AS)
N_AS = N_AS0 / math.exp(-r_AS * T_EU_AS)
# Migration rates during the various epochs.
m_AF_B = 25e-5
m_AF_EU = 3e-5
m_AF_AS = 1.9e-5
m_EU_AS = 9.6e-5
# Population IDs correspond to their indexes in the popupulation
# configuration array. Therefore, we have 0=YRI, 1=CEU and 2=CHB
# initially.
population_configurations = [
msprime.PopulationConfiguration(
sample_size=0, initial_size=N_AF),
msprime.PopulationConfiguration(
sample_size=1, initial_size=N_EU, growth_rate=r_EU),
msprime.PopulationConfiguration(
sample_size=1, initial_size=N_AS, growth_rate=r_AS)
]
migration_matrix = [
[ 0, m_AF_EU, m_AF_AS],
[m_AF_EU, 0, m_EU_AS],
[m_AF_AS, m_EU_AS, 0],
]
demographic_events = [
# CEU and CHB merge into B with rate changes at T_EU_AS
msprime.MassMigration(
time=T_EU_AS, source=2, destination=1, proportion=1.0),
msprime.MigrationRateChange(time=T_EU_AS, rate=0),
msprime.MigrationRateChange(
time=T_EU_AS, rate=m_AF_B, matrix_index=(0, 1)),
msprime.MigrationRateChange(
time=T_EU_AS, rate=m_AF_B, matrix_index=(1, 0)),
msprime.PopulationParametersChange(
time=T_EU_AS, initial_size=N_B, growth_rate=0, population_id=1),
# Population B merges into YRI at T_B
msprime.MassMigration(
time=T_B, source=1, destination=0, proportion=1.0),
# Size changes to N_A at T_AF
msprime.PopulationParametersChange(
time=T_AF, initial_size=N_A, population_id=0)
]
# Use the demography debugger to print out the demographic history
# that we have just described.
dp = msprime.DemographyDebugger(
Ne=N_A,
population_configurations=population_configurations,
migration_matrix=migration_matrix,
demographic_events=demographic_events)
dp.print_history()
def variable_recomb_example():
infile = "../hapmap/genetic_map_GRCh37_chr22.txt"
# Read in the recombination map using the read_hapmap method,
recomb_map = msprime.RecombinationMap.read_hapmap(infile)
# Now we get the positions and rates from the recombination
# map and plot these using 500 bins.
positions = np.array(recomb_map.get_positions()[1:])
rates = np.array(recomb_map.get_rates()[1:])
num_bins = 500
v, bin_edges, _ = scipy.stats.binned_statistic(
positions, rates, bins=num_bins)
x = bin_edges[:-1][np.logical_not(np.isnan(v))]
y = v[np.logical_not(np.isnan(v))]
fig, ax1 = pyplot.subplots(figsize=(16, 6))
ax1.plot(x, y, color="blue")
ax1.set_ylabel("Recombination rate")
ax1.set_xlabel("Chromosome position")
# Now we run the simulation for this map. We assume Ne=10^4
# and have a sample of 100 individuals
tree_sequence = msprime.simulate(
sample_size=100,
Ne=10**4,
recombination_map=recomb_map)
# Now plot the density of breakpoints along the chromosome
breakpoints = np.array(list(tree_sequence.breakpoints()))
ax2 = ax1.twinx()
v, bin_edges = np.histogram(breakpoints, num_bins, density=True)
ax2.plot(bin_edges[:-1], v, color="green")
ax2.set_ylabel("Breakpoint density")
ax2.set_xlim(1.5e7, 5.3e7)
fig.savefig("_static/hapmap_chr22.svg")
if __name__ == "__main__":
# single_locus_example()
# multi_locus_example()
# mutations_example()
# segregating_sites_example(10, 5, 100000)
# migration_example()
# out_of_africa()
variable_recomb_example()