add example of "add more mutations to SLiM sim"

[petrelharp]

As @bhaller reports, doing

initialize() {
	initializeTreeSeq();
	initializeMutationRate(1e-12);
	initializeMutationType("m1", 0.5, "f", 0.0);
	initializeMutationType("m2", 0.5, "f", 0.01);
	initializeGenomicElementType("g1", m2, 1.0);
	initializeGenomicElement(g1, 0, 1e8-1);
	initializeRecombinationRate(1e-8);
}
1 early() {
	sim.addSubpop("p1", 500);
}
5000 late() {
	sim.treeSeqOutput("overlay.trees");
}

then

import msprime, tskit

ts = tskit.load("overlay.trees").simplify()
mutated = msprime.sim_mutations(ts, rate=1e-7, random_seed=1, keep=True)
mutated.dump("final_overlaid.trees")

results in

Traceback (most recent call last):
 File "/Users/bhaller/Desktop/ts_overlay.py", line 4, in <module>
  mutated = msprime.sim_mutations(ts, rate=1e-7, random_seed=1, keep=True)
 File "/Users/bhaller/Library/Python/3.9/lib/python/site-packages/msprime/mutations.py", line 1692, in sim_mutations
  _msprime.sim_mutations(
msprime._msprime.LibraryError: Existing allele(s) incompatible with mutation model alphabet.

which is confusing. The solution is to add

ts = pyslim.convert_alleles(pyslim.generate_nucleotides(ts))

but this is not well-explained anywhere. We should add an example to the docs of this.

[bhaller]

Huh, interesting. Why was this problem not reproducible for me? I ran it over and over, after I got the initial error, and was never able to reproduce the problem. Also, is this a new wrinkle? I've had an example like this for quite a long time, and it has never been a problem until now.

[bhaller]

Huh, interesting. Why was this problem not reproducible for me? I ran it over and over, after I got the initial error, and was never able to reproduce the problem. Also, is this a new wrinkle? I've had an example like this for quite a long time, and it has never been a problem until now.

Your discussion on Slack perhaps answered these questions for me; the problem is rare, for this model, because it happens only when a JK mutation occurs on an edge that a SLiM mutation is above. So, hmm. Kinda hoses the simple example I thought I had for this workshop worksheet. :-<

[bhaller]

OK, so. It would indeed be great to have this discussed in the pyslim docs, and once that happens, I'll link to that from the SLiM Workshop's worksheet on tree sequence recording, which has an example that occasionally hits this problem.

Clemens Weiss mentioned on Slack:

passing msprime.SLiMMutationModel as the mutation model to msprime.sim_mutations should do the trick, as shown e.g. here: https://tskit.dev/pyslim/docs/latest/tutorial.html#adding-neutral-mutations-to-a-slim-simulation

That would be a natural place to add a mention of ts = pyslim.convert_alleles(pyslim.generate_nucleotides(ts)) as an alternative fix, probably; so, no need for a whole new section, perhaps, but rather just a note about this alternative approach. Right now that section says:

We’ll add these using the msprime.SLiMMutationModel, so that the file can be read back into SLiM, but any of the other mutation models in msprime could be used.

I guess that comes the the caveat that ts = pyslim.convert_alleles(pyslim.generate_nucleotides(ts)) might be needed for those other mutational models to work; so that would be a natural place to add this, maybe?

Anyhow, I've got my brain wrapped around this now, and will decide what to do for the example in the workshop. :->