Msprime: Allow binary ancestral and derived states

[petrelharp]

In SLiM we are recording derived states as a sequence of integers, packed in binary-wise. This creates problems in the python interface, for instance:

tabs = msprime.TableCollection(sequence_length=1.0)
tabs.nodes.add_row()
tabs.nodes.add_row(time=1.0)
tabs.edges.add_row(left=0.0, right=1.0, parent=1, child=0)
tabs.sites.add_row(position=0.5, ancestral_state='')
tabs.mutations.add_row(site=0, node=0, derived_state=b'\x96\x00\x00\x00\x00\x00\x00\x00')

print(tabs)

In this case, printing throws an error, but the problem occurs not only with printing.

One possible solution is to always use the errors="surrogateescape" option to the utf-8 codec, as described here and in PEP-0383, which allows lossless round-tripping with unprintable characters. I've tried this out over in #543. However, this doesn't seem great, because (a) it doesn't work, and I am having trouble finding where the encode is that still throws an error; and (b) it is new to python3.

The subject of encodings is still pretty opaque to me. Do you have any suggestions, @jeromekelleher?

[jeromekelleher]

That's a bit tricky @petrelharp, I hadn't thought people were going to pack binary data into these... being able to put arbitrary binary data in there basically makes it impossible to provide nice printing options. This is going to be text data after all in the vast majority of cases/

Would it work at all if you put in (say) a comma separated list of the integers you want to represent as strings? Tskit is never going to intrepret these values anyway, so I don't see any significant performance hit implied in this. Conversion could be done at the end during the output phase in SLiM so the cost of generating the strings would be negligible.

[petrelharp]

We could do that. However, I'd naively thought that if we we're printing anything, then we wouldn't run into the issue. Here's an example of the error occurring without anything being printed:

# as above, then
ts = tabs.tree_sequence()

for m in ts.mutations():
    pass

What do you think, @bhaller?

[jeromekelleher]

We could fix the underlying problem with the encodings @petrelharp, but l'm not sure it's the right thing to do from the end-user's perspective. I think it's a lot more user-friendly for them to be able to print out sites and mutations as expected without having to specify obscure encoding options to various methods.

Sorry for my ignorance here, but what do the integers mean? How big are they, and what is the expected length of the lists? Are they interpretable outside of the context of a SLiM invocation?

[petrelharp]

Every new mutation in SLiM gets an ID, which is an int64_t; since these can "stack", i.e., a genome can carry more than one of these at once at a site, we record the derived_state to a new mutation as being the current list of SLiM mutation IDs, and stick the list of structs recording information about these (selection coef, etc) in the metadata. To think concretely about it, imagine you've got two sets of mutation types, one is ACGT, the other is methylation; these naturally stack.

So yeah: this does not lend itself to pretty-printing. An end-user might want to do something with it, but it'll be on pyslim to provide an interface to that.

Here's another option: part of the pre-processing that pyslim does is to take the list of derived states, store them somewhere else, and then replace them with some text (probably, 0, 1, 2, ...). That way, printing haplotypes and the like will work just fine; and the information will not be lost.

One slightly wierd thing about that plan is that we don't currently have a python method to load in just the tables; we only have msprime.load() that makes the whole tree sequence. In this case it'd be nice to load in just the tables; then modify them; and finally make a tree sequence.

[jeromekelleher]

Here's another option: part of the pre-processing that pyslim does is to take the list of derived states, store them somewhere else, and then replace them with some text (probably, 0, 1, 2, ...). That way, printing haplotypes and the like will work just fine; and the information will not be lost.

Sounds like a good plan!

One slightly wierd thing about that plan is that we don't currently have a python method to load in just the tables; we only have msprime.load() that makes the whole tree sequence. In this case it'd be nice to load in just the tables; then modify them; and finally make a tree sequence.

Adding TableCollection.load would be straightforward now --- it's really just a case of some simple plumbing between the Python layer and the C API. Definitely one to add soon.

It's not strictly necessary to do this though, is it? We can call msprime.load() and get the tables OK so long as we don't try to do any processing, right?

[petrelharp]

It's not strictly necessary to do this though, is it? We can call msprime.load() and get the tables OK so long as we don't try to do any processing, right?

Nope; it's fine as-is. But I think we should do that plumbing sooner or later. Maybe I'll have a go at it.

[bhaller]

Hey all, just catching up on this thread. So it seems like we have various options:

1. SLiM could write out ASCII data for derived states, converting in from binary on the way out and converting it from ASCII on the way in

2. pyslim could do that same conversion on the Python side, allowing the file format to stay as it is

3. The derived state data could be kept as binary, and tskit could be fixed not to raise when that is the case. Printing would still not work, unless some way of providing tskit with a custom print-formatter from pyslim, or some such mechanism, were devised.

I see why printing raises, but why do other tskit operations raise? What do non-printing methods do with the derived-state data that depends upon it being ASCII?

It seems like it would be possible for tskit to make printing safe too; it could scan the derived state for non-ASCII bytes, and if any are present, fall back upon just printing raw byte values or some such. This would not make the output useful, but it would prevent the crash. But pyslim could just provide a substitute print method, to be used when the derived state data is known to be from SLiM, no? If that is workable, that would seem to be the simplest and perhaps best solution. And maybe tskit could somehow know to automatically use the pyslim printer, or at least be told to do so by pyslim, so that the right printing would be done automatically when a SLiM file was used (at least when pyslim was explicitly brought into the process of reading it)?

Any of 1, 2, 3 seem workable; this seems like it's mostly a matter of policy. I would prefer not to do (1) because we're very close to shipping SLiM 3 and I don't want to risk a major bug in that release. To me, (3) seems like the right course of action: proclaim that derived states can contain binary data, as a matter of policy, and fix things up so that that works without crashing, and then figure out a solution for the printing problem. If we decide that derived state data must be ASCII, that will not only require a conversion to ASCII when SLiM is saving to disk; it will also require a conversion whenever SLiM builds a tree sequence (during coalescence-checking, during runtime cross-checking). That will introduce overhead that seems to me to be unnecessary. In fact, if the policy is that derived-state data must be ASCII, then to be compliant with policy SLiM would even need to convert to ASCII before every simplification, no? I really don't think we want to go down that road; that would be an enormous upheaval to SLiM's code, days before shipping the introductory release of this new feature. Instead, I think derived state data should be regarded the same way metadata is: client-specific binary.

[petrelharp]

Hm, good points. A couple of things:

1. the problem with non-ASCII derived states is purely python, not in C; so you don't have to worry about it in cross-checking.

2. I think the correct way to do (3) would be something like: look through the ancestral/derived states, and if they are >1 character, assign them digits (or letters or something), and then print those instead of the actual state, having some way to look up what they actually are.

3. But, that's exactly how we'd want to do (2); only doing it in "general" would be substantially harder because it'd be hard to anticipate what people would want to do in general. So, that's why I'm voting for option (2).

[bhaller]

the problem with non-ASCII derived states is purely python, not in C; so you don't have to worry about it in cross-checking.

OK, as long as we all agree that, as a matter of policy, one can have binary data in derived states in C code and on disk, just not in Python. If that's agreed, then that's fine.

I think the correct way to do (3) would be something like: look through the ancestral/derived states, and if they are >1 character, assign them digits (or letters or something), and then print those instead of the actual state, having some way to look up what they actually are.

Aha. This sounds quite reasonable.

But, that's exactly how we'd want to do (2); only doing it in "general" would be substantially harder because it'd be hard to anticipate what people would want to do in general. So, that's why I'm voting for option (2).

Hmm. OK, so it seems like we're choosing between (2) and (3). I don't have a strong feeling, between the two of them; (3) strikes me as a bit more "pure" of a solution, but it sounds like you're arguing that (2) would be easier to implement/use? The only thing that gives me pause about (2), really, is that it seems to imply that usage of pyslim for reading in SLiM .trees files is mandatory. Reading them without using pyslim would mean that binary derived-state data would be present, which would be a violation of policy and would cause unpredictable errors (as it does now, I guess). Using pyslim would be required, to patch that up on the way in and out. Is that where we want to be? Seems a bit unfortunate; (3) seems better to me because it seems much simpler to be able to say "you don't have to use pyslim just to read in a SLiM-generated .trees file, you only need to use it if you want to print mutations in a comprehensible format, or if you want to write out SLiM-compliant files". But if you guys feel strongly that (2) is better, I don't oppose it; you understand these Python issues far better than I, and you're the ones writing the code, so I'm happy to defer to you.

[petrelharp]

OK, as long as we all agree that, as a matter of policy, one can have binary data in derived states in C code and on disk, just not in Python. If that's agreed, then that's fine.

I agree with this.

The only thing that gives me pause about (2), really, is that it seems to imply that usage of pyslim for reading in SLiM .trees files is mandatory.

Good point, but I don't think it's a problem to require pyslim, as you'll need that to get some of the information out of the .trees file anyhow. But, I do think it would be good to make it so python does not error on the tree sequences as long as the offending binary states are not printed; this would make pyslim not mandatory. (I don't know what's involved in doing that, though.)

(3) strikes me as a bit more "pure" of a solution, but it sounds like you're arguing that (2) would be easier to implement/use?

Currently we have no use case for this re-labeling of derived states outside of SLiM, so it's moot whether to put the functionality in pyslim or tskit. If it were in tskit, I do think it should be as an optional method, not as something that happens by default (you don't want programs automatically relabeling alleles!). So, how about I'll implement it in pyslim, and if it looks general/useful enough, we can move it up into tskit at some point?

[bhaller]

This all sounds good to me. So (2), but with the element from (3) that tskit will not error out on binary derived states except in printing, if that is possible (or perhaps printing can be protected from causing an error too, as I proposed above, but since the output would be gibberish anyway maybe who cares). OK. Let me know if there's anything I can do to help, but since this is all python stuff I don't imagine there is...

[jeromekelleher]

OK, as long as we all agree that, as a matter of policy, one can have binary data in derived states in C code and on disk, just not in Python. If that's agreed, then that's fine.

I'm becoming more convinced that I don't agree with this, unfortunately! My reasoning does revolve around Python, but I think that it applies more generally too as more languages handle unicode better.

In Python 3+, the str type represents unicode text data which is utf-8 (a superset of ASCII) encoded by default, and supports all sort of other unicode encodings. String data is unicode at a fundamental level. The bytes type on the other hand represents arbitrary binary data. You can store unicode or ASCII data in there just fine, but if you ever want to turn it into a string (i.e., for printing), you must decode it under some encoding.

In tskit ancestral and derived states are unicode strings. There is a very specific and concrete reason for this choice. Users will want to do things like:

for site in ts.sites():
    if site.ancestral_state == "A":
        # do something

If ancestral and derived states were bytes rather than strings we would have to write this as:

for site in ts.sites():
    if site.ancestral_state == b"A":
        # do something

Note the bytes literal b"A" now. So fine, this will work, but it will be a source of terrible bugs because users will naturally think that writing site.ancestral_state == "A" will work, and will not expect to have to put in the cryptic b prefix. Worse, it's very, very easy to forget the b prefix leading to silent and pernicious bugs in user code.

So, this is the fundamental reason that I chose to represent ancestral and derived states as unicode strings, not as arbitrary bytes: library users would --- justifiably! --- hate us. It's worth repeating here that one of the key applications of tskit is going to be working with real data, which will overwhelmingly consist of simple character data for the ancestral and derived states. Users of the library in this context will absolutely want to perform text manipulation and comparisons on the states, and they should be able to do this in as natural and un-error-prone way as possible.

The reason that operations other than printing fail when non-unicode data is present in the states is that every time we create a Site or Mutation object, we also create a unicode str object representing the ancestral/derived state. This happens in all sorts of places, and so the changes required to hack around having non-unicode data in there would be deep and fundamental.

So, after thinking about this some more, I think that we should require that ancestral and derived states must be unicode in all conforming tree sequences. Otherwise we'll end up in all sorts of messes. So, under this interpretation there is only one solution to the present problem: SLiM must convert the binary data to unicode/utf-8/ascii before outputting the final tree sequence.

I don't want to be too prescriptive about this, but I do feel strongly about site.ancestral_state == "A" always working as expected, and I don't see the value in the very tricky mucking about behind the scenes that would be needed to support this and the possibility of having non-unicode data in the states. Surely all of unicode is enough!

[bhaller]

Well, this is what I was worried about. SLiM won't only need to convert the binary data to ASCII before outputting the final tree sequence, though; it will need to do it before making a tree sequence in all situations, such as when doing a runtime cross-check, or when checking for coalescence. This is assuming that it is permissible to have arbitrary binary data for the derived states when it is recording tables; if that is not permitted, either, then it would need to convert to ASCII every time a new mutation is recorded. That was quite a large performance hit when we were doing it, though, so I hope we don't need to go back to that.

I do wonder, though, if it is necessary to be quite so strict about this. The issues raised all seem to have to do (a) with Python, and (b) with other users of the library. Well, in Python I think we have agreed that pyslim will always insulate tskit from the binary data; if a SLiM .trees file is being read in, pyslim will convert the binary data to ASCII before a tree sequence is ever made, as I understand it. And other users of the library are, of course, free to always place unicode data into the derived state, so their usage pattern would, perhaps, never be violated. The cases that you are concerned about would simply never come up, it seems to me. (If we are worried about users opening a SLiM .trees file without using pyslim, I suppose we could devise a safeguard specifically to prevent that!)

The question, I think, is whether it is necessary for any of the core code inside tskit, in C, to know whether the derived state data is binary or unicode. Is that, in itself, a problem? Or do the concerns really just revolve around the representation when in Python, which can be managed by pyslim? You say:

The reason that operations other than printing fail when non-unicode data is present in the states is that every time we create a Site or Mutation object, we also create a unicode str object representing the ancestral/derived state. This happens in all sorts of places, and so the changes required to hack around having non-unicode data in there would be deep and fundamental.

It isn't clear to me whether you're speaking here of Python, or C, or both, however.

Anyway, if this is how it needs to be, it's not the end of the world. It will add extra runtime overhead in SLiM whenever cross-checking or coalescence-checking are enabled, but those features are not enabled by default. Unfortunately it will mean another change to the file format, and all of my saved .trees files will be invalid once again; but better to sort that out before we ship SLiM 3. Let's settle this ASAP, and I'll patch up SLiM if that's what needs doing.

[jeromekelleher]

I do wonder, though, if it is necessary to be quite so strict about this. The issues raised all seem to have to do (a) with Python, and (b) with other users of the library. Well, in Python I think we have agreed that pyslim will always insulate tskit from the binary data; if a SLiM .trees file is being read in, pyslim will convert the binary data to ASCII before a tree sequence is ever made, as I understand it. And other users of the library are, of course, free to always place unicode data into the derived state, so their usage pattern would, perhaps, never be violated. The cases that you are concerned about would simply never come up, it seems to me. (If we are worried about users opening a SLiM .trees file without using pyslim, I suppose we could devise a safeguard specifically to prevent that!)

We shouldn't need pyslim to intervene when loading tree sequences --- it should be possible to load any tree sequence that is stored on file and perform basic operations on it. OK, you might not be able to interpret the metadata, but you should be able to perform general operations without it (computing stats, etc). This is basic interoperability and it's vital: suppose I had a whole pile of tree sequence files that have come from various different source (msprime, tsinfer, SLiM, ...), and I needed to upgrade the file format. I want to be able to run msp upgrade on all of these. Having to figure out which one is a SLiM tree sequence so I can massage it then load it with the correct library is a huge loss in generality. (OK in this example you probably could just msp upgrade all the files and it'd work fine, but you get the point.)

It isn't clear to me whether you're speaking here of Python, or C, or both, however.

The C library doesn't know anything about unicode at the moment, so arbitrary bytes can be stored.

Anyway, if this is how it needs to be, it's not the end of the world. It will add extra runtime overhead in SLiM whenever cross-checking or coalescence-checking are enabled, but those features are not enabled by default. Unfortunately it will mean another change to the file format, and all of my saved .trees files will be invalid once again; but better to sort that out before we ship SLiM 3. Let's settle this ASAP, and I'll patch up SLiM if that's what needs doing.

I don't think it's necessarily a run-time cost for SLiM. Until you actually output the tree sequence file, you can store whatever you want in the ancestral and derived states. I don't think we'll ever check for unicode correctness at that level, as it'd be too complicated and expensive. You only need to do these conversions at the end, just before writing to a file. I can't see this being a significant overhead.

So, to be clear, my suggestion is that we require that all tree sequence files must have ancestral and derived states encoded as UTF-8.

[bhaller]

So, to be clear, my suggestion is that we require that all tree sequence files must have ancestral and derived states encoded as UTF-8.

OK. So in SLiM, both when using tables and when building tree sequences from tables, the derived state can be binary? That does improve the situation, yes. I misunderstood you before when you wrote:

I think that we should require that ancestral and derived states must be unicode in all conforming tree sequences.

I took that to mean in C as well as in Python and on disk. But as clarified now, I'm happy to agree; this seems entirely reasonable. I agree with your interoperability point, it does seem best for .trees files to be fundamentally readable by tskit, always, without needing an auxiliary library to translate things. @petrelharp, are you on board?

[jeromekelleher]

OK. So in SLiM, both when using tables and when building tree sequences from tables, the derived state can be binary? That does improve the situation, yes. I misunderstood you before when you wrote:

Yes, when working with the C API basically anything goes for these states, and I think this is the correct approach (C is inherently dangerous, after all). Once the tree sequence files are written to disk the states must be UTF-8 encoded if they are to be interoperable. We'll need to document this somewhere...

As we're going through this it's made me realise what's important and not, so this has been a really useful discussion; thanks!

[petrelharp]

@petrelharp, are you on board?

I'm on board, and it seems for the best: the proposal is to have SLiM, when it writes out tree sequence files, convert the ancestral states of the Site table and the derived states of the Mutation table to ASCII, right?

[bhaller]

OK; so let it be written, so let it be done. I'll fix SLiM today.

[bhaller]

OK, I have committed a fix; SLiM now writes (and reads) ASCII derived state data. I'm running my full test suite on it now, but I'm pretty sure it's working properly. Peter, a reminder: I'm working on the master branch of the Messer Lab SLiM repository now. I'm not sure what branch you're working with now, but you might need to merge changes over.

@petrelharp wrote above:

convert the ancestral states of the Site table and the derived states of the Mutation table to ASCII

I'm not doing anything with the ancestral states in the site table, because they're empty in SLiM anyway, right? None of our existing ASCII-fication code does anything with the site table. Let me know if that needs to change!

[petrelharp]

I'm not doing anything with the ancestral states in the site table, because they're empty in SLiM anyway, right? None of our existing ASCII-fication code does anything with the site table. Let me know if that needs to change!

I forget if simplification might cause a derived state to be moved to the ancestral state (e.g., if there's a mutation on the root branch) - isn't that an option, maybe? If so, you'd need to check the ancestral states also.

I'll close this issue now: we've decided WONTFIX NOTABUG.

[bhaller]

I forget if simplification might cause a derived state to be moved to the ancestral state (e.g., if there's a mutation on the root branch) - isn't that an option, maybe? If so, you'd need to check the ancestral states also.

Well, in SLiMSim::TreeSequenceDataToAscii() there's a comment, written by you :->, that says:

/*** Notes: ancestral states are always zero-length, so we don't need to Ascii-ify Site Table ***/

And we don't do anything with ancestral states anywhere in SLiM as far as I know. If simplification does that, that would be a problem; but since my test suite all runs without errors, I suspect it does not.