diff --git a/examples/examples.py b/examples/examples.py index e1e1501..d496c38 100644 --- a/examples/examples.py +++ b/examples/examples.py @@ -45,6 +45,7 @@ help="name of file to output selected locus information", default="sel_loci.txt") parser.add_argument("--treefile","-t", type=str, dest="treefile", help="name of output file for trees (default: not output)",default=None) +parser.add_argument("--seed", "-d", dest="seed", type=int, help="random seed") args = parser.parse_args() @@ -56,6 +57,8 @@ import simuOpt simuOpt.setOptions(alleleType='mutant') import simuPOP as sim +sim.setRNG(seed=args.seed) +random.seed(args.seed) def fileopt(fname,opts): '''Return the file referred to by fname, open with options opts; @@ -135,14 +138,17 @@ def __call__(self, loc, alleles): id_tagger = sim.IdTagger() id_tagger.apply(pop) first_gen = pop.indInfo("ind_id") -init_ts = msprime.simulate(2*len(first_gen), - length=max(locus_position)) +length = max(locus_position) + +# Since we want to have a finite site model, we force the recombination map +# to have exactly `length` loci with a fixed recombination rate between them. +rcmb_map = msprime.RecombinationMap.uniform_map(length, args.recomb_rate, length) +init_ts = msprime.simulate(2*len(first_gen), recombination_map=rcmb_map) haploid_labels = [(k,p) for k in first_gen for p in (0,1)] node_ids = {x:j for x, j in zip(haploid_labels, init_ts.samples())} rc = RecombCollector(ts=init_ts, node_ids=node_ids, locus_position=locus_position) - # initially, population is monogenic init_geno=[sim.InitGenotype(freq=1.0)] diff --git a/examples/proc.py b/examples/proc.py new file mode 100644 index 0000000..76e39bb --- /dev/null +++ b/examples/proc.py @@ -0,0 +1,8 @@ +# coding: utf-8 +import msprime +ts = msprime.load('neutral_ts') +tsc = msprime.TreeStatCalculator(ts) +print('mean pairwise tmrca matrix among all 6 haploid samples for 10 windows on genome') +for t in tsc.mean_pairwise_tmrca(windows=[float(i) for i in range(10)], + sample_sets=[[0,1,2,3,4,5], [0,1,2,3,4,5]]): + print(t) diff --git a/examples/test.sh b/examples/test.sh new file mode 100644 index 0000000..5e636f3 --- /dev/null +++ b/examples/test.sh @@ -0,0 +1,6 @@ +echo "trees differ across genome with recombination" +python ./examples.py -d 111 -T 1 -N 3 -r 1e-1 -L 10 -a 0.0001 -b 0.0001 -k 3 -t neutral_ts -g /dev/null > /dev/null && python proc.py +echo " " +echo "trees conserved across genome without recombination" +python ./examples.py -d 111 -T 1 -N 3 -r 0 -L 10 -a 0.0001 -b 0.0001 -k 3 -t neutral_ts -g /dev/null > /dev/null && python proc.py + diff --git a/examples/tree_comparison.txt b/examples/tree_comparison.txt new file mode 100644 index 0000000..6b55dca --- /dev/null +++ b/examples/tree_comparison.txt @@ -0,0 +1,132 @@ +simuPOP Version 1.1.8.3 : Copyright (c) 2004-2016 Bo Peng +Revision 4553 (Feb 11 2017) for Python 3.5.3 (64bit, 1thread) +Random Number Generator is set to mt19937 with random seed 0x62d96efe2e082f18. +This is the standard mutant allele version with 256 maximum allelic states. +For more information, please visit http://simupop.sourceforge.net, +or email simupop-list@lists.sourceforge.net (subscription required). +Options: +Namespace(chrom_length=10, gamma_alpha=0.0001, gamma_beta=0.0001, generations=10000, logfile='-', neut_mut_rate=1e-07, nsamples=5000, nselloci=1, outfile=None, popsize=2, recomb_rate=0.1, sel_mut_rate=1e-07, selloci_file='sel_loci.txt', simplify_interval=500, treefile='neutral_ts') +10:09:07 11/01/17 PDT +---------- + + Init ts +############################################################ +# Nodes # +############################################################ +id flags population time +0 1 0 0.00000000000000 +1 1 0 0.00000000000000 +2 1 0 0.00000000000000 +3 1 0 0.00000000000000 +4 0 0 0.24138545192840 +5 0 0 0.36768962082101 +6 0 0 0.80818994340547 +7 0 0 1.10295116208317 +8 0 0 1.21146637134289 +9 0 0 1.55374021332996 +10 0 0 1.90995962338469 +11 0 0 3.09020402491767 +12 0 0 4.12376873296128 +############################################################ +# Edges # +############################################################ +id left right parent child +0 0.00000000 9.00000000 4 2 +1 0.00000000 9.00000000 4 3 +2 2.00000000 3.00000000 5 1 +3 2.00000000 3.00000000 5 4 +4 0.00000000 2.00000000 6 0 +5 0.00000000 2.00000000 6 4 +6 2.00000000 4.00000000 7 0 +7 3.00000000 4.00000000 7 4 +8 2.00000000 3.00000000 7 5 +9 8.00000000 9.00000000 8 0 +10 0.00000000 2.00000000 8 1 +11 8.00000000 9.00000000 8 4 +12 0.00000000 2.00000000 8 6 +13 7.00000000 8.00000000 9 1 +14 7.00000000 8.00000000 9 4 +15 4.00000000 7.00000000 10 1 +16 4.00000000 7.00000000 10 4 +17 4.00000000 5.00000000 11 0 +18 3.00000000 4.00000000 11 1 +19 3.00000000 4.00000000 11 7 +20 4.00000000 5.00000000 11 10 +21 5.00000000 8.00000000 12 0 +22 8.00000000 9.00000000 12 1 +23 8.00000000 9.00000000 12 8 +24 7.00000000 8.00000000 12 9 +25 5.00000000 7.00000000 12 10 +############################################################ +# Sites # +############################################################ +id position ancestral_state +############################################################ +# Mutations # +############################################################ +id site node derived_state parent +############################################################ +# Migrations # +############################################################ +id left right node source dest time + + In recomb collector +############################################################ +# Nodes # +############################################################ +id flags population time +0 1 0 0.00000000000000 +1 1 0 0.00000000000000 +2 1 0 0.00000000000000 +3 1 0 0.00000000000000 +4 0 0 0.24138545192840 +5 0 0 0.36768962082101 +6 0 0 0.80818994340547 +7 0 0 1.10295116208317 +8 0 0 1.21146637134289 +9 0 0 1.55374021332996 +10 0 0 1.90995962338469 +11 0 0 3.09020402491767 +12 0 0 4.12376873296128 +############################################################ +# Edges # +############################################################ +id left right parent child +0 0.00000000 9.00000000 4 2 +1 0.00000000 9.00000000 4 3 +2 2.00000000 3.00000000 5 1 +3 2.00000000 3.00000000 5 4 +4 0.00000000 2.00000000 6 0 +5 0.00000000 2.00000000 6 4 +6 2.00000000 4.00000000 7 0 +7 3.00000000 4.00000000 7 4 +8 2.00000000 3.00000000 7 5 +9 8.00000000 9.00000000 8 0 +10 0.00000000 2.00000000 8 1 +11 8.00000000 9.00000000 8 4 +12 0.00000000 2.00000000 8 6 +13 7.00000000 8.00000000 9 1 +14 7.00000000 8.00000000 9 4 +15 4.00000000 7.00000000 10 1 +16 4.00000000 7.00000000 10 4 +17 4.00000000 5.00000000 11 0 +18 3.00000000 4.00000000 11 1 +19 3.00000000 4.00000000 11 7 +20 4.00000000 5.00000000 11 10 +21 5.00000000 8.00000000 12 0 +22 8.00000000 9.00000000 12 1 +23 8.00000000 9.00000000 12 8 +24 7.00000000 8.00000000 12 9 +25 5.00000000 7.00000000 12 10 +############################################################ +# Sites # +############################################################ +id position ancestral_state +############################################################ +# Mutations # +############################################################ +id site node derived_state parent +############################################################ +# Migrations # +############################################################ +id left right node source dest time