Save a mask in the tree sequence

[hyanwong]

In 2a7341c two sorts of mask were implemented for VCF output. In some cases, we might want to save a default mask within the tree sequence itself, for example one that masks off the telomeres, or (in the case of the tree sequence of, say, the q arm of a chromosome) a large region which is there simply in order to get the positions aligned correctly. By default, these regions could also be ignored by the stats framework (see #346 for a related discussion). @gtsambos noticed this problem with the unified genealogy tree sequences recently.

Would it be possible to store a mask somehow at the top level of the tree sequence. I assume this isn't metadata, if it is used by functions such as write_VCF or the stats stuff.

[gtsambos]

It sounds like there are two different types of information that potentially rely on this operation:

Summary information about lengths
For instance, sequence_length. To me it would be most obvious if this referred to the total length of the unmasked bits of sequence. And if this isn't possible, there should be some easy way for users to figure out what the length of the unmasked sequence is because this will often be needed as a normalisation factor in downstream calculations.

Information about positions
To me, it would make most sense if this corresponded to the actual position along the chromosome (as I think is the case already, if I've understood Yan correctly)

So for instance, say you have a chromosome of length 100 with the p arm going from 0-45, the centromere going from 45-55, and the q arm from 55-100. Then if tsq is the tree sequence holding the q arm, I think it would make most sense if tsq.sequence_length was 45, with tsq.first() having a left coordinate of 55 and tsq.last() having a right coordinate of 99.

[hyanwong]

Hmm, I think changing ts.sequence_length so that it did not correspond to the last position along the sequence would cause a lot of breakage? It seems like we want a different property, e.g. unmasked_span (or something like that) which we would use in e.g. the stats framework? In cases where there was no mask, that would be equal to sequence_length, of course.

[hyanwong]

Oh, and I do think that sequence_length is not a great name, and a better one would be sequence_end because (a) it conflates the position and the extent of the sequence (as @gtsambos points out above) and (b) we are trying to use "span" rather than "length" to avoid confusion with branch lengths in trees (which measure time).

Nevertheless, sequence_length is so embedded in tskit that I don't think we are likely to want to rename it. At least, not without unanimous positive support.

[benjeffery]

There's no changing the meaning or name of sequence_length at this point, for sure!

[benjeffery]

I'm not sure about implicit use of a stored mask, it wouldn't be obvious to the reader of some code that one was being used. It would also be a breaking change to all methods that you modified to use one, as the current expectation is that all sites will be used.

Once we have mask usage in the stats framework it may be that storing a mask would be useful, at least in top-level metadata.

[hyanwong]

For the stats framework, it may be enough to ignore any regions in which there is no topology and no sites? Perhaps ts.nonmissing_span (or ideally a better name) could return the span for which this is not the case?

[petrelharp]

I think that we probably should try out some unofficial ways of storing masks first (so, in our workflows store masks in top-level metadata, etc), and add options for masks to things like the statistics, and after use cases are more clear, then think about how/whether to make a new slot in TreeSequence for them.