Using pyslim to join msprime and slim simulation for samples early in time

[freshspaceoctopus]

Hello, I would like to start a slim simulation from a msprime simualted tree with some samples that were sampled early in time (sample time is 400 unlike other samples). For those specific samples, I would like to simulate a WF in SLiM while the other samples that were sampled at time 0 in msprime can wait. The only way I could think of in my head is to 1. simplify that specific samples -> change time 400 to 0 -> import into slim after applying pyslim annotation.

Many thanks!

[bhaller]

Hi @freshspaceoctopus! I'm not sure you need to change the time to 0. SLiM can start a simulation at whatever tick you like. Just do your load of the saved .trees file in tick 400 (or maybe 401?). But maybe I'm not getting a forward-time vs. backward-time wrinkle here that makes it more difficult than that. @petrelharp will have more intelligent things to say about that than I do. :->

[petrelharp]

I more or less agree, @bhaller, but you're missing one detail: when you load in the .trees file, it will set the tick counter to whatever the tick counter was when you saved it. (I think this is what you meant, but might be confusing.) And, there's an argument tick (and cycle) to pyslim.annotate that lets you set this value.

I think this will take some messing around to get it to work, however.

[bhaller]

Ah, that is true. Well, you can change the tick/cycle counters in SLiM to be whatever you want; they are writeable properties. So maybe that's the easy approach? Changing the counters after you load into SLiM? But I agree, it might take some messing around to find the best solution – and to be quite certain it is working correctly.

[freshspaceoctopus]

Thank you @bhaller and @petrelharp for your reply,
Sorry for the ambiguous post, I was in a hurry when I was writing the original post, so it was too brief.
Let me explain my problem in depth.

For example, I have a tree file simulated in msprime like this:

import tskit, msprime

# Demography
demography = msprime.Demography()
demography.add_population( name = "p0", initial_size = 10000)
demography.add_population( name = "p1", initial_size = 10000, default_sampling_time = 500)
demography.add_population( name = "p0_p1", initial_size = 10000)
demography.add_population_split(time = 1000, ancestral = "p0_p1", derived = ["p0", "p1"])

# Sampling
sample_list = []
sample_list.append(msprime.SampleSet(20, population = "p0")
sample_list.append(msprime.SampleSet(10000, population = "p1")

# Recombination rate
r_chrom = 1e-6
r_break = math.log(2)
chrom_positions =  np.arange(0, 10e6 + 1, 1e6)
map_positions = []
rates = []
for i in range(len(chrom_positions)):
        map_positions.append(chrom_positions[i])
        if i > 0:
            rates.append(r_chrom)
        if i > 0 and i+1 < len(chrom_positions):  # Add the +1 element for positions 1, 2, and 3
            map_positions.append(chrom_positions[i] + 1)
            rates.append(r_break)
rate_map = msprime.RateMap(position=map_positions, rate=rates)



# Simulation
ts = msprime.sim_ancestry(sample_list,
                              demography = demography,
                              model = [
                                msprime.DiscreteTimeWrightFisher(duration = 10),
                                msprime.StandardCoalescent(duration = 490),
                                msprime.DiscreteTimeWrightFisher(duration = 10),
                                msprime.StandardCoalescent()],
                              random_seed = 1,
                              recombination_rate = rate_map,
                              ploidy = 2)

The demography itself is a simplified version, but what I want to do is create a burn-in for slim using msprime. I will provide neutral mutations, then use a subset of these neutral mutations that p0 and p1 share to simulate quantative trait selection for 500 generations for subpopulation "p1".
So, the sampling time of p1 is earlier than p0 in the msprime simulation: 500 generations back in the past from the present. I aim to fill the 500 generations with SLiM simulation.

One approach I used was to modify the time value of the tree sequence so that I only need to load in "p1" when simulating the 500 generations of forwards-in time using SLiM, and I can just use the command pyslim.annotate with tick being 1, but this felt a little bit hacky, since I have to modify the original tree table. Since SLiM starts ticks at 1, I was wondering if I can provide negative values such as -500 so that SLiM can start only on the simulation of p1, and then when it reachs 0 or 1, it can continue on with both p0 and p1.

This is off topic, but the reason I specified the msprime model as above was because I wanted to simulate discrete wright fisher for 10 generations for each subpopulation, then follow up with the Hudson coalescent, but didn't know how to approach it with subpopulations with default_sampling_time > 0.

Another solution that comes in mind is to set merging time of p0 and p1 to 500 generations in msprime, sprinkle neutral mutations for SLiM's selection, simulate another 500 generations in msprime for p0, simulate 500 generations in SLiM for p1, and use tskit.union to merge all three trees. Would this be the standard way of tackling this issue?

Many thanks in advance!

(edit) Edited minor errors

[freshspaceoctopus]

So, I kinda overcome the above issue by pulling the original tree 500 generations early (so that there is a negative value of time in the tree_sequence), use pyslim's simplify to extract only the p1 individuals, run SLiM, then load it back into python and use tskit.decapitate, then use tskit.union to merge it back to the original msprime simulated tree. Maybe I am overcomplicating my approach, but at least it seems to work...

[petrelharp]

Hi again, @freshspaceoctopus! Thanks for the additional details; they confirm that we were interpreting your original question correctly.

So, I think the tldr; is still that there isn't a standard way to do what you're asking - that is, it's a reasonable thing to do, but I don't know of anyone having done it. Making it work is going to involve (a) a bit of messing around, and (b) some verification that it's doing what you want.

And - what you describe in the last comment -- simplify -> SLiM -> union -- sounds like a relatively elegant way to deal with it. Nice job working that out. To make sure, I'd do some good checking that things are as you expect: for instance, that MRCAs between your final p1 and p2 genomes are all in the time range you'd expect.

[freshspaceoctopus]

Thanks, @petrelharp for the confirmation!
After checking MRCA and connection of nodes, it seems to be working quite well!
Thank you for your time, and thank everyone for tskit, msprime, SLiM and pyslim: these are truely amazing!