Tskit: simplify removes an ancestor of a sample (!?!?)

[petrelharp]

A report on the SLiM mailing list by Megan Smith has turned out out to be a bug in simplify (!?!?!) as far as I can tell.

With this SLiM script as sim.slim (it's pretty quick):

initialize() {
    initializeTreeSeq();
    initializeMutationRate(7.5e-7);
    initializeMutationType("m1", 0.25, "g", -0.0133, 0.35);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 9999);
    initializeRecombinationRate(5e-6);
}
1 {
    sim.addSubpop("p1", 1250);
}
12500 { 
    sim.addSubpopSplit("p2", 1250, p1); 
    sim.treeSeqRememberIndividuals(sim.subpopulations.individuals);
}
13750 late() {
    sim.treeSeqOutput("sim.trees");
}

running this python script:

import os
import pyslim, tskit
import numpy as np

os.system("slim -s 1809121419536 sim.slim > sim.log")

ts = pyslim.load("sim.trees")

keep_nodes = [9994, 9995, 7698, 7699, 8020, 8021, 8314, 8315, 8810, 8811, 9638, 9639, 9294, 9295, 9654, 9655, 7856, 7857, 9944, 9945, 12298, 12299, 11576, 11577, 12200, 12201, 11466, 11467, 12100, 12101, 12350, 12351, 10866, 10867, 10244, 10245, 12342, 12343, 10274, 10275]

sts = ts.simplify(keep_nodes)

nid = 14
on = ts.node(keep_nodes[nid])
n = sts.node(nid)
assert(on.metadata.slim_id == n.metadata.slim_id)

ot = ts.at(1652)
t = sts.at(1652)
print(f"The parent at position 1652 of this node:")
print(on)
print(f"in the original tree is:")
print(ts.node(ot.parent(on.id)))
print(f"but in the simplified tree, that node is:")
print(n)
print(f"with parent {t.parent(14)}.")
print(t.draw(format='unicode'))

tells us that indeed, the 14th sample in the simplified tree sequence has no parent for a small chunk of the chromosome, when it did in the original tree sequence. That's bad!

The parent at position 1652 of this node:
{'id': 9654, 'time': 0.0, 'population': 1, 'individual': 1077, 'metadata': NodeMetadata(slim_id=37504654, is_null=0, genome_type=0), 'flags': 1}
in the original tree is:
{'id': 12839, 'time': 1.0, 'population': 1, 'individual': -1, 'metadata': NodeMetadata(slim_id=37498685, is_null=0, genome_type=0), 'flags': 0}
but in the simplified tree, that node is:
{'id': 14, 'time': 0.0, 'population': 0, 'individual': 7, 'metadata': NodeMetadata(slim_id=37504654, is_null=0, genome_type=0), 'flags': 1}
with parent -1.

The bug only showed up due to missing genotypes in the output, and Megan reports it happened around 5% of the time - not a one-off, apparently (but it does depend on which nodes are simplified for this tree sequence).

This is not a pyslim bug; the same thing happens if we tskit.load( ) (I've used pyslim to show the exta information).

I currently have no idea what's going on.

[jeromekelleher]

Crud. Are you tracking this one down ~Alistair~ Peter, or will I look into it? Might be worth pushing the ts through the python version, if it's not too big?

[petrelharp]

I assume you meant "Peter" instead of "Alistair"? If you could look into it that'd be wonderful; I'm supposed to be cutting a page out of a grant application atm. The ts is only 22249 nodes and 34953 edges - should be do-able.

[jeromekelleher]

Hah, yes, Peter. My multitasking is pretty buggy! I'll take a look ASAP.

[molpopgen]

Is it possible that this is some truly odd topology? I hope I did this right, but I attempted to trace all parental nodes from the node in question to root, and then ask if the paths to root for all other nodes crossed it. I added this to the end of the script:

u = on.id
path_to_root = []
while u != tskit.NULL:
    # if u != on.id:
    path_to_root.append(u)
    u = ot.parent(u)

print(path_to_root)

YES=0
for k in keep_nodes:
    if k != on.id:
        u = k
        while u != tskit.NULL:
            # print(u)
            if u in path_to_root:
                YES += 1
            u = ot.parent(u)

print(f"Times we crossed the path = {YES}")

The answer:

Times we crossed the path = 0

So, that suggests tossing in a keep_unary=True, and the (abridged) output changes to:

The parent at position 1652 of this node:
{'id': 9654, 'time': 0.0, 'population': 1, 'individual': 1077, 'metadata': NodeMetadata(slim_id=37504654, is_null=0, genome_type=0), 'flags': 1}
in the original tree is:
{'id': 12839, 'time': 1.0, 'population': 1, 'individual': -1, 'metadata': NodeMetadata(slim_id=37498685, is_null=0, genome_type=0), 'flags': 0}
but in the simplified tree, that node is:
{'id': 14, 'time': 0.0, 'population': 0, 'individual': 7, 'metadata': NodeMetadata(slim_id=37504654, is_null=0, genome_type=0), 'flags': 1}
with parent 46.

Perhaps we need an annual/quarterly award for who can generate the longest unary node path?

[molpopgen]

Oh, if you draw an SVG of the original tree in question, there are 2 roots, so it hasn't be recapitated. I'd post it, but it is big.

[molpopgen]

Seems that recapitation does fix the issue. Doing a recapitate with the same recombination rate as the slim script, removing my change to keep unary nodes gives me:

The parent at position 1652 of this node:
{'id': 9654, 'time': 0.0, 'population': 1, 'individual': 1077, 'metadata': NodeMetadata(slim_id=37504654, is_null=0, genome_type=0), 'flags': 1}
in the original tree is:
{'id': 12839, 'time': 1.0, 'population': 1, 'individual': -1, 'metadata': NodeMetadata(slim_id=37498685, is_null=0, genome_type=0), 'flags': 0}
but in the simplified tree, that node is:
{'id': 14, 'time': 0.0, 'population': 0, 'individual': 7, 'metadata': NodeMetadata(slim_id=37504654, is_null=0, genome_type=0), 'flags': 1}
with parent 848.

[jeromekelleher]

I'm with @molpopgen here - sure looks like there's just two roots in the first tree, and that sample has a nice long unary path back to one of them:

ot = ts.at(1652)
print("old tree roots = ", ot.roots)
for node in keep_nodes:
    u = node
    while ot.parent(u) != -1:
        u = ot.parent(u)
    print(node, "->", u)

t = sts.at(1652)
print("new tree roots", t.roots)
u = n.id
while u != -1:
    print(u)
    u = t.parent(u)

gives

old tree roots =  [19, 2075]
9994 -> 19
9995 -> 19
7698 -> 19
7699 -> 19
8020 -> 19
8021 -> 19
8314 -> 19
8315 -> 19
8810 -> 19
8811 -> 19
9638 -> 19
9639 -> 19
9294 -> 19
9295 -> 19
9654 -> 2075
9655 -> 19
7856 -> 19
7857 -> 19
9944 -> 19
9945 -> 19
12298 -> 19
12299 -> 19
11576 -> 19
11577 -> 19
12200 -> 19
12201 -> 19
11466 -> 19
11467 -> 19
12100 -> 19
12101 -> 19
12350 -> 19
12351 -> 19
10866 -> 19
10867 -> 19
10244 -> 19
10245 -> 19
12342 -> 19
12343 -> 19
10274 -> 19
10275 -> 19
new tree roots [14, 707]
14

[molpopgen]

Glad this seems to be it. It's not even 10AM and I need a beer after reading this.

[jeromekelleher]

I've unlabelled this as a "bug", but I'll let you close it @petrelharp.

[petrelharp]

Hah, of course! Thanks a bunch, @molpopgen.

[jeromekelleher]

Glad this seems to be it. It's not even 10AM and I need a beer after reading this.

Haha, I hear ya! One thing I could totally do without is a bug in simplify right now!