Less memory-intensive method for getting a single haplotype, including ancestral haplotypes

[hyanwong]

A number of use-cases have come up recently where people want to get a single ancestral haplotype out of a tree sequence. The quick-and-dirty way to do this, as used by @saurabhbelsare , is to change the flag on the ancestral node to mark it as a sample, and then use the haplotypes iterator. This (presumably) will report the unknown sections as "missing data". I wonder if we want a .haplotype() method that gets a single haplotype (possibly an ancestral one) from a tree sequence.

Alternatively, and this is maybe better, we could wrap the method suggested above by modifying the .haplotypes() iterator to (optionally) take a list of nodes: haplotypes(nodes=None, *, ...). If nodes is None, we default to the sample nodes, as currently, but if nodes is not None it can be a list of any nodes, including non-sample ones. In this case, we take a copy of the tables, mark the specified nodes as samples (and all other nodes as non-samples), simplify(), and output the haplotypes. This would also be a way to allow people to output a few sample haplotypes without slurping the whole variant matrix into memory (at the moment we advise that if only a few haplotypes are needed, the user should simplify first, but we don't actually do the work for them).

At the moment all the parameters to ts.haplotypes() are keyword only, so this won't break API compatibility.

[jeromekelleher]

No reason we can't do this, I think. Already possible with the variants method, AFAIK.

[hyanwong]

Thansk @jeromekelleher. I think the variants() method only allows sample nodes to be specified. I could either adjust the variants() method to allow non-sample nodes too, in which case haplotypes becomes a simple wrapper, or implement the "change node flags + simplify" part within haplotypes() in which case variants would still only allow sample nodes to be used. Preferences?

[jeromekelleher]

method only allows sample nodes to be specified.

Pretty sure it allows non-sample "samples"

[hyanwong]

Ah, perhaps the docs need changing then. I didn't actually check the C code. Will do.

(Edit - this is is the docs, sorry - I was confused because of the parameter naming)

[hyanwong]

variants() does allow any node (not just samples), but requires isolated_as_missing=False otherwise we get

_tskit.LibraryError: Cannot generate genotypes for non-samples when isolated nodes are considered as missing

So it's not so helpful for creating ancestral haplotypes (in which I assume we want to represent the missing chucks with a missing_data_character rather than the ancestral state)

[jeromekelleher]

Yes, that's true. Not sure there's a simple way to do what you want then, detecting missingness requires the node in question to be a sample.

[hyanwong]

It works if, in the haplotypes() function, we detect when there are non-sample nodes with isolated_as_missing=True and simply create an identical TS where those nodes are marked as samples. I have this in a PR - I guess I can submit it, although it seems a little bit of a hack to me.

[petrelharp]

Is there a fundamental reason that missingness only applies to samples? It seems perfectly consistent that isolated non-samples be treated as missing.

[hyanwong]

Is there a fundamental reason that missingness only applies to samples? It seems perfectly consistent that isolated non-samples be treated as missing.

Yes, I was wondering that too. I hadn't got very far in my thinking, though.

[hyanwong]

We should probably close this in favour of #1207, unless we want to finish any discussion on missingness in non-samples.

[jeromekelleher]

Is there a fundamental reason that missingness only applies to samples? It seems perfectly consistent that isolated non-samples be treated as missing.

I think we'd need to do some refactoring in the tree generation code to allow "samples" to be set dynamically for this to work - it's not something I'm enthusiastic about taking on though, as that code is subtle and it would need a lot of testing.

[hyanwong]

Also see #3298 for an alternative unimplemented approach.