Running into a new error due to mismatch between Expected root time and Observed root times

[elissasoroj]

Hello!

Per usual, I have been running simulations incorporating alternation of generations, which stores the haploid and diploid individuals in separate populations. For post-sim analysis these individuals need to end up back in the same population. I used to use the following code to move all the individuals into one population ("p0") then I would use simplify() with keep_input_roots=True to remove the now-empty population:

#copy the tables
tables = tree_sequence.dump_tables()

#change the population to p0 for everyone
nnodes = tables.nodes.time.size
tables.nodes.population = np.zeros(nnodes, dtype=np.int32)

# drop nodes that are not connected to anything.
nodes_in_edge_table = list(set(tables.edges.parent).union(tables.edges.child))

# remove the empty population nodes by using simplify
tables.simplify(
    samples=nodes_in_edge_table,
    keep_input_roots=True,
    filter_individuals=True,
    filter_populations=True,
    filter_sites=False
)

# turn it back into a treesequence
new_ts = tables.tree_sequence()

Recently this functionality has broken, I suspect due to a pyslim update. I've tried changing various setting in the simplify() function, as well as a different method for selecting the samples (I tried keeping them all), to no avail. It seems like some of the pertinent nodes are being removed because some of the root times are no longer the ancestral root time (see error).

Error, for reference:

ValueError: Not all roots of the provided tree sequence are at the time expected by recapitate(). This could happen if you've simplified in python before recapitating (fix: don't simplify first). If could also happen in other situations, e.g., you added new individuals without parents in SLiM during the course of the simulation with sim.addSubPop(), in which case you will probably need to recapitate with msprime.sim_ancestry(initial_state=ts, ...). (Expected root time: 201; Observed root times: [34.0, 201.0, 43.0, 46.0, 17.0, 27.0])

If I don't run the simplify() function then the recapitation works just fine. My impression is that using keep_input_roots=True should prevent this issue. What am I missing?

The empty population may not seem like a big deal, but it ends up being a nuisance for more complex situations where demography is involved. A workaround may be editing the tables manually to remove the empty population, but alas I couldn't figure out how to do that.

model as .txt - script_ex.txt

Thank you for maintaining and updating pyslim!
~Elissa

[petrelharp]

Hello! Thanks for the easy-to-run code.
Let's see - if you really just want to only remove that entry in the population table, you can just do like this:

tables = tree_sequence.dump_tables()
#change the population to p0 for everyone
nnodes = tables.nodes.num_rows
tables.nodes.population = np.zeros(nnodes, dtype=np.int32)
tables.populations.clear()
tables.populations.append(tree_sequence.population(0))
# turn it back into a treesequence
new_ts = tables.tree_sequence()

pyslim.recapitate(new_ts, ancestral_Ne=10)

That works (for me).

Why are you getting that error? Something weird is happening with your simplify-but-mark-everyone-as-samples-first call. I haven't dug into why exactly this happens, but I think you probably don't want to do it, anyhow? I have confirmed that if you don't do that then it also works:

tables = tree_sequence.dump_tables()
nnodes = tables.nodes.num_rows
tables.nodes.population = np.zeros(nnodes, dtype=np.int32)
tables.populations.clear()
tables.populations.append(tree_sequence.population(0))

tables.simplify(
    keep_input_roots=True,
    filter_individuals=True,
    filter_populations=True,
    filter_sites=False
)

new_ts = tables.tree_sequence()
pyslim.recapitate(new_ts, ancestral_Ne=10)

This also works (for me).

If you do want to simplify something else away, explain what it is you're wanting to do?

FYI, the only thing (I think) that has changed recently in pyslim is the consistency check that you're hitting - so, it might well be saving you from something nonsensical happening.

[petrelharp]

Okay, I looked into this, and this happened because there were isolated, non-sample nodes in your original tree sequence (i.e., nodes with no parent and no child, but that were not marked as samples); this is fine, but then your python script marked those as samples. Since 'roots' are defined to be highest-up places accessible from samples, these weren't roots before, but then became roots. Here's the code (which also checks that in fact - as one would expect if simplifying while marking everything as a sample - all the nodes and all the trees are the same):

edge_nodes = np.unique(np.concatenate([tree_sequence.edges_parent, tree_sequence.edges_child]))
ts = tree_sequence.simplify(
    edge_nodes,
    keep_input_roots=True,
)
pyslim.recapitate(ts, ancestral_Ne=10)

np.all(ts.tables.nodes.metadata_vector("slim_id") == tree_sequence.tables.nodes.metadata_vector("slim_id"))
ts.tables.edges == tree_sequence.tables.edges

bad_roots = []
for t in ts.trees():
  for r in t.roots:
     if ts.node(r).time < 201:
        bad_roots.append((t.index, r))


t0 = tree_sequence.first()
t1 = ts.first()
for k, n in bad_roots:
    t0.seek_index(k)
    t1.seek_index(k)
    assert(t0.parent(n) == -1)
    assert(t1.parent(n) == -1)
    assert(t0.left_child(n) == -1)
    assert(t1.left_child(n) == -1)
    print("isolated node:", n, "flags:", tree_sequence.node(n).flags, ts.node(n).flags)

As for why those are in the file, I'm not sure - you do have simplification turned off; if not, that would definitely have removed these already. Ah, I think I see - I'm guessing it has something to do with the haploid phases and null genomes? I see things like this:

            child = p1.addRecombinant(individual.genome1, NULL, NULL,
                         individual.genome2, NULL, NULL, parent1 = individual, parent2 = individual);

and

            child = p0.addRecombinant(individual.genome1, NULL, NULL,
                          NULL, NULL, NULL, parent1 = individual); //only 1 parent recorded

... so, maybe those are the null genomes? I don't remember how those work under the hood. It might be worth double-checking that this isn't a symptom of something going wrong. (However, the fact that the ancestries of all the genomes alive at the end of the simulation are fine, and go back to the start of the simulation suggests everything is fine?)

Edit: they aren't null genomes, according to SLiM:

>>> ts.node(n)
Node(id=1746, flags=1, time=26.0, population=1, individual=-1, metadata={'slim_id': 34795129, 'is_null': False, 'genome_type': 0})

[bhaller]

Huh. Null genomes seemed like it made so much sense. Hopefully we're setting that is_null flag correctly everywhere.

[elissasoroj]

I do think there was an error in my models that didn't mark two parents for every diploid child. I'm not sure what happens under the hood (I think we discussed it and thought it didn't matter), but I have changed the models to mark parent1 and parent2, even if it is the same individual (but to be clear, there were some lines that didn't mark both parents correctly when they were different individuals).
This may or may not fix the issue. If you want me to do some testing let me know.

[bhaller]

Hmm. I don't think the parent1/parent2 parameters to addRecombinant() ought to matter to this. The tree sequence records the 0/1/2 genomes that the new gamete inherits from, regardless of those parameters. Those parameters are just for pedigree tracking, I think. But I would like to know where these stray genomes are coming from, if they are not null genomes! How can we track that down?

[elissasoroj]

I can do some digging!

[bhaller]

Programming rule of thumb: The more consistency checks the better. :->