tskit-dev / pyslim

Tools for dealing with tree sequences coming to and from SLiM.
MIT License
27 stars 23 forks source link

Cannot replicate the example in extracting-information-about-selected-mutations #204

Closed xin-huang closed 2 years ago

xin-huang commented 2 years ago

Hi,

I tried to run the example in https://pyslim.readthedocs.io/en/latest/tutorial.html#extracting-information-about-selected-mutations

I simulated the data with the codes from https://pyslim.readthedocs.io/en/latest/tutorial.html#extracting-information-about-selected-mutations

initialize()
{
   initializeSLiMModelType("WF");
   initializeTreeSeq();
   initializeMutationRate(1e-6);
   initializeMutationType("m1", 0.5, "e", -0.1);
   initializeMutationType("m2", 0.5, "e", 0.5);
   initializeGenomicElementType("g1", c(m1, m2), c(0.9, 0.1));
   initializeGenomicElement(g1, 0, 1e6-1);
   initializeRecombinationRate(1e-8);
}

1 early() {
   sim.addSubpop("p1", 1000);
}

1000 {
   sim.treeSeqOutput("selection.trees");
}

and set the seed to 23 as suggested in the tutorial.

Here is the output from SLiM 3.6

// Initial random seed:
23

// RunInitializeCallbacks():
initializeSLiMModelType(modelType = 'WF');
initializeTreeSeq();
initializeMutationRate(1e-06);
initializeMutationType(1, 0.5, "e", -0.1);
initializeMutationType(2, 0.5, "e", 0.5);
initializeGenomicElementType(1, c(m1, m2), c(0.9, 0.1));
initializeGenomicElement(g1, 0, 999999);
initializeRecombinationRate(1e-08);

// Starting run at generation <start>:
1 

Then I loaded the simulated data

ts = pyslim.load("selection.trees")
print(ts.num_mutations)
# 5961
print(ts.num_sites)
# 5941

But my result is ts.num_mutations = 6044 and ts.num_sites = 6020

If I print ts.mutation(0), my result is

Mutation(id=0, site=0, node=3010, derived_state='1653896', parent=-1, metadata={'mutation_list': [{'mutation_type': 2, 'selection_coeff': 1.5596729516983032, 'subpopulation': 1, 'slim_time': 827, 'nucleotide': -1}]}, time=172.0)

While in the tutorial

m = ts.mutation(0)
print(m)
# {'id': 0,
#  'site': 0,
#  'node': 4425,
#  'time': 1.0,
#  'derived_state': '1997240',
#  'parent': -1,
#  'metadata': {
#     'mutation_list': [
#           { 'mutation_type': 2,
#             'selection_coeff': 1.4618088006973267,
#             'subpopulation': 1,
#             'slim_time': 992,
#             'nucleotide': -1
#           }
#     ]
#  }
# }
xin-huang commented 2 years ago

The results are the same as in https://tskit.dev/pyslim/docs/latest/tutorial.html#extracting-information-about-selected-mutations

petrelharp commented 2 years ago

Ok, glad it's all good!