Existing alleles incompatible with mutation model alphabet

[mshpak76]

After recapitating a slim-generated tree using pyslim, I use sim_mutations in msprime to map neutral mutations on the recapitated tree. The slim tree is generating binary 0/1 alleles at the site(s) undergoing selective sweep(s), and in most instances the following sequence of functions generates a plausible mutated tree sequence:

ts = tskit.load("ThreePopFix.trees")
ts_new = pyslim.SlimTreeSequence(ts)
rts = pyslim.recapitate(ts_new, recombination_rate = 1.14e-8, gene_conversion_rate = 5.9e-8, gene_conversion_tract_length = 518, ancestral_Ne=20000, random_seed = ancestry_seed)
mts = msprime.sim_mutations(rts, rate = 1.22e-8, keep=True, random_seed = mutation_seed, model=msprime.BinaryMutationModel())

However, in some instances msprime.sim_mutations() returns the following error:

_msprime.sim_mutations(
msprime._msprime.LibraryError: Existing allele(s) incompatible with mutation model alphabet.

Strangely, in many cases executing the sim_mutations() code a second (or third) time resolves the problem and generates output. I can't discern any pattern for when it does and doesn't work, particularly since re-running it often creates output, and the fact that the pyslim output format is the same in all cases.

[petrelharp]

Hi, @mshpak76! This is because mutations generated by SLiM do not have 0/1 alleles - SLiM needs to store a lot more information about the mutations. This is covered, for instance, here. The reason you get the error only sometimes is that sometimes when you add mutations, no mutations land on previously-mutated sites, and so sim_mutations never notices that there's existing alleles it doesn't know how to make sense of.

So, probably what you want to do is to add mutations using the SLiM Mutation Model. However, you should think carefully about what's going on here - if you add more mutations to an already-mutated site under selection, that means that there were neutral mutations at that same site as well that didn't affect the process of selection. This might confuse your downstream analysis.

[mshpak76]

Is there a way to do this that prevents pyslim from mapping mutations onto the sites that have been mutated in slim? In my case this should be straightforward enough, I'm using slim to model a selective sweep so that the only mutations introduced into the ts by slim are the site(s) under selection.

Alternatively, is there a setting in Slim to force a compatible allele format output in the ts?

[petrelharp]

Yes - use a RateMap for sim_mutations() that excludes the selected site.

As for the second question: you could use a nucleotide model in SLiM and then use pyslim.convert_alleles( ). You'd not be able to load the result back in to SLiM, though.

[mshpak76]

Thanks, this is very helpful!