Annotating neutral coalescent with multiple subpops

[dylandebaun]

Hello!
I'm hoping to load in a neutral simulation from msprime into SliM. I am not sure how to scale up the vignette to work with multiple populations.

I get the following error in SliM:
sim.readFromPopulationFile("ts_annotated.init.trees");
ERROR (Species::__TabulateSubpopulationsFromTreeSequence): node is_vacant unexpected; this file cannot be read.

After running this code to produce the tree sequence.

import msprime, pyslim

seed = 1
Ne = 10000
time_units = 2 * Ne

demography = msprime.Demography()
for pop in ["n0", "n1", "n2", "n3", "n4", "n5", "n6"]:
    demography.add_population(name=pop, initial_size=Ne)

scale = -2.0
demography.add_population_split(time=(2.5 + scale) * time_units, derived=["n1", "n2"], ancestral="n4")
demography.add_population_split(time=(3.5 + scale) * time_units, derived=["n4", "n0"], ancestral="n5")
demography.add_population_split(time=(7 + scale) * time_units, derived=["n5", "n3"], ancestral="n6")

demography.sort_events()

breaks = [0, 3e6, 6e6, 9e6]
recomb_map = msprime.RateMap(position=breaks, rate=[5e-8, 0.5e-8, 5e-8])

# Simulate ancestry
ots = msprime.sim_ancestry(
    samples={"n0": 100, "n1": 100, "n2": 100, "n3": 100},  
    demography=demography,
    recombination_rate=recomb_map,
    random_seed=seed
)

ots = pyslim.annotate(ots, model_type="WF",tick=1, stage="late")

ots=pyslim.remove_vacant(ots)

tables = ots.tables
ts_metadata = tables.metadata
ts_metadata["SLiM"]["model_type"] = "WF"
tables.metadata = ts_metadata
ots = tables.tree_sequence()

ots.dump("ts_annotated.init.trees")

[petrelharp]

I believe the problem is with this line:

ots=pyslim.remove_vacant(ots)

"removing vacant nodes" (really, removing the sample flag from them) is something you want to do for analysis in python; however, to load the tree sequence into SLiM you need to include them.

So, I suspect this will work if you delete that line?

I also don't think you should have this code chunk:

tables = ots.tables
ts_metadata = tables.metadata
ts_metadata["SLiM"]["model_type"] = "WF"
tables.metadata = ts_metadata
ots = tables.tree_sequence()

The annotate line already specifies model_type="WF", so this code chunk should do nothing at all.

I'm curious if there's somewhere in the docs that suggested you needed to do those steps?

Please go ahead and close if this solves the issue.

[dylandebaun]

Hi Peter,
I obtained that code chunk from this part of the vignette and wasn't sure if it was required for all simulations: https://tskit.dev/pyslim/docs/stable/vignette_coalescent_diversity.html#load-into-slim

Removing that code chunk and ots=pyslim.remove_vacant(ots) produces the same error:

import msprime, pyslim

seed = 1
Ne = 10000
time_units = 2 * Ne

demography = msprime.Demography()
for pop in ["n0", "n1", "n2", "n3", "n4", "n5", "n6"]:
    demography.add_population(name=pop, initial_size=Ne)

scale = -2.0
demography.add_population_split(time=(2.5 + scale) * time_units, derived=["n1", "n2"], ancestral="n4")
demography.add_population_split(time=(3.5 + scale) * time_units, derived=["n4", "n0"], ancestral="n5")
demography.add_population_split(time=(7 + scale) * time_units, derived=["n5", "n3"], ancestral="n6")
demography.sort_events()

breaks = [0, 3e6, 6e6, 9e6]
recomb_map = msprime.RateMap(position=breaks, rate=[5e-8, 0.5e-8, 5e-8])
# Simulate ancestry
ots = msprime.sim_ancestry(
    samples={"n0": 100, "n1": 100, "n2": 100, "n3": 100},  
    demography=demography,
    recombination_rate=recomb_map,
    random_seed=seed
)

ots = pyslim.annotate(ots,model_type="WF",tick=1,stage="late")
# ots=pyslim.remove_vacant(ots)
# tables = ots.tables
# ts_metadata = tables.metadata
# ts_metadata["SLiM"]["model_type"] = "WF"
# tables.metadata = ts_metadata
# ots = tables.tree_sequence()

ots.dump("ts_annotated.init.trees")

In Slim:

initialize()
{
    defineConstant("L", 9000000);
    initializeSLiMModelType("WF");
    initializeTreeSeq(timeUnit="generations");
    initializeMutationRate(c(0.03e-8, 0.003e-8, 0.03e-8),
                           c(3000000, 6000000, 9000000) - 1);
    initializeMutationType("m1", 0.5, "f", 0.0);
    initializeMutationType("m2", 0.5, "e", 0.01);
    initializeGenomicElementType("g1", m2, 1.0);
    initializeGenomicElement(g1, 0, L-1);
    initializeRecombinationRate(1e-8);
}


1 late() {
	// Import msprime tree sequence
	sim.readFromPopulationFile("ts_annotated.init.trees");

ERROR (Species::__TabulateSubpopulationsFromTreeSequence): node is_vacant unexpected; this file cannot be read.

[petrelharp]

Oh darn; I was sure that was the problem.

But - I've just run the code above and everything works fine? Are you sure you're using the newly produced tree sequence? With SLiM v5.1?

Thanks very much for the working, quick-to-run example, btw.

[dylandebaun]

Ah I see what happened, I needed to upgrade tskit 0.6.0-py313hde99e4e_0 --> 0.6.4-py313hc577518_1
Thanks so much for sanity checking this for me!