Tskit: Edge case issue with simplify and missing data

[hyanwong]

I've just hit an edge case which is preventing me round-tripping some missing data examples. If we have an extreme edge of the genome in which only a single sample has non-missing data, then this can be represented by a tree at that point with only a single branch, connecting that sample to the root. However, if we run simplify() on such a tree sequence, the edge is removed (as it only contains unary nodes). That leaves the sample as an "isolated node", and hence the missing data code in https://github.com/tskit-dev/tskit/pull/272/ flags it up as a case where the genotype should be set to -1, even though in this case, we do have information to properly encode the genotype.

I'm wondering if this is a issue with the missing data code, or the simplify() code? For example, in simplify() it might be considered reasonable not to drop unary nodes from a sample if they connect that sample to the root? But I'm not sure how the root would be identified in this case.

Ping @jeromekelleher and @petrelharp as they are the simplifying and missing data gurus :)

[jeromekelleher]

Hmm. This is a good corner case I hadn't thought of. We have two samples in the tree, one is missing data and the other isn't. How do we encode this? This is tricky --- the current isolated roots criteria is good because it's very cheap to do. More complicated rules are going to get expensive and will probably end up with worse corner cases.

[petrelharp]

Well, this is arguing for encoding missingness explicitly rather than through topology. But it's a very minor corner case, so no need to ditch the whole concept. I think that pushes the problem over to simplify. Maybe we just need to code this corner case into simplify? Hmph.

[hyanwong]

FWIW this feels like a corner case for simplify() to me. In this case the unary connection between the sample and the root is meaningful: it definitely reflects a genealogical relationship that we might want to keep, and not simplify() away.

The fact that it then provides a way to correctly report sequence data under the "isolated nodes == missing" criterion merely reflects the underlying genealogical truth.

[petrelharp]

Hm: I think that simplify is definately doing the right thing, as originally defined. That edge isn't reflecting a genealogical relationship between the samples, which is how we've defined things.

[jeromekelleher]

I agree with Peter here. I guess you can specify the keep_unary option in simplify in this case (which we will want to do in tsinfer anyway; otherwise we split up a bunch of edges), and the problem goes away. It is a pretty obscure corner case - if data is missing for n - 1 of your samples, do you really know much about the n'th?

[petrelharp]

if data is missing for n - 1 of your samples, do you really know much about the n'th?

Oh, good point. If only one is nonmissing, then probably that one should be missing also.

[hyanwong]

Fair points. I'm (possibly wrongly) thinking of the root as somehow special, a bit like a sample. It does seem a little different from a "normal" unary node anyhow.

Oh, good point. If only one is nonmissing, then probably that one should be missing also.

Some of this comes down again to whether we want to encode the sample data itself as a tree sequence, in which case we want to be able to recover what we put in. I do imagine that there might be cases where only one of the samples has data (e.g. if this is a subset of a larger number of samples).

Either way, for the moment I'll just pass the "keep unary" option, I guess. I don't know if this is an argument for changing keep_unary=T|F to something like keep="unary"|"root"|None

[petrelharp]

Thinking a bit more about this with @andrewkern: it might be helpful to think about two kinds of missing data:

1. missing topology
2. missing genotypes

Disconnecting a sample from the tree indicates the first thing. The proposal to, in the future, designate a "missing" derived state could be used to indicate the second thing. If we're always inferring topology from genotypes, then as @hyanwong points out, there's always going to be an arbitrary decision whether to call a bit of the topology missing or not, but that's OK and analogous to the decision to call a genotype missing or not.

So, is simplify doing the right thing? It seems like it is creating a chunk of missing topology that was not missing before. But, simplify only cares about relationships between samples, and in the input tree sequence there were no relationships between samples, so it correctly marks all topology as missing. Note that if the root was a sample, then this would not happen.

[hyanwong]

simplify only cares about relationships between samples

I think this is what confuses me (a little). I expected simplify to care about relationships between samples and between samples and the root (but not care about other non-sample nodes). But that's probably just me.

[petrelharp]

I think this is what confuses me (a little). I expected simplify to care about relationships between samples and between samples and the root (but not care about other non-sample nodes). But that's probably just me.

Well, "the root" is dependent on the samples: if you add a sample that coalesces above the current MRCA of all the previous samples, you have to move the root up. Conversely, removing samples might move the MRCA (i.e., the root of the tree that relates them) down. Here, I'm equating "root" with "MRCA". You maybe have a different notion?

[hyanwong]

I have the same notion of the root as you, and it is indeed a function of the samples themselves. I suspect I just haven't thought through the meaning of simplify as thoroughly as you have.

[jeromekelleher]

I think we've reached the conclusion that this is an odd corner case, but (a) it's a pretty extreme one that's unlikely to happen; and (b) it will automatically get handled in the only case that we do know it happens. So, I think we can close this issue?

[petrelharp]

I agree, although it was useful to think through.

[hyanwong]

Actually, I think it probably will happen a lot in the example I'm thinking about: where we treat each partial fragment from a sequencer as a separate sample. In this case, with a small number of samples (low coverage) it's pretty likely that there will be places on the left or right where there is only coverage from one sample.

But I can deal with this using simplify(keep_unary=True). I'd prefer a simplify(keep_root=True) option, but it seems a little excessive to code this up for this particular corner case. If there are other reasons why such a flag to simplify would be useful, I'd be interested to hear them, though.

[andrewkern]

Actually, I think it probably will happen a lot in the example I'm thinking about: where we treat each partial fragment from a sequencer as a separate sample. In this case, with a small number of samples (low coverage) it's pretty likely that there will be places on the left or right where there is only coverage from one sample.

strikes me that simplify is not what you would want to do in this case and instead you want to treat these data as missing genotypes as @petrelharp mentions above

[jeromekelleher]

But I can deal with this using simplify(keep_unary=True). I'd prefer a simplify(keep_root=True) option, but it seems a little excessive to code this up for this particular corner case. If there are other reasons why such a flag to simplify would be useful, I'd be interested to hear them, though.

But we want to have keep_unary=True as the default in tsinfer anyway, so this really isn't an issue. We'll compress the data much better with keep_unary=True, I can explain why in person when I get back.

[hyanwong]

OK - I didn't think about this! I guess therefore for my purposes this is moot, and this issue can be closed?