SLiM needs a place to put mutation metadata

[bhaller]

Since the dawn of tree-sequence recording in SLiM, we've been putting metadata about mutations (selection coefficient, mutation type id, origin subpopulation, origin tick, nucleotide) into the derived state information. That has never been great. One reason is that since the derived state has to be ASCII, we have to convert all this metadata to ASCII on the way out, and back to binary on the way back in, which is slow and messy. Another reason is that since SLiM can "stack" mutations, creating a new derived state that includes the previous mutation at the site as well as a new mutation "stacked" on top of it, we end up re-emitting all of this information whenever that happens (which is rare in many models, but commonplace in some models). So we end up with multiple copies of what might be duplicate information about a given mutation – or might even not be duplicate information, since the state of the mutation might have been changed in the meantime. A third reason why this is not great is that the final state of the mutation is not captured; if the user changes the mutation type of a mutation and then writes out, that state change is (I'm 99% sure?) not captured, since no new derived state gets recorded. In all, we've got a bloated, inefficient, confusing mess.

But we've been limping along with this. The straw that breaks the camel's back is that I am now shifting SLiM to support multiple traits. A given mutation can have pleiotropic effects on any/all traits in the model, which means that it needs to have separate selection and dominance coefficients for each trait in the model. So the potential amount of metadata per mutation, and the complexity of it, is about to jump upward sharply. Even in present-day SLiM people have been doing models with something like six or eight traits (with zero support in tree-sequence recording for that, and close to zero support in SLiM for it too!); once doing this is properly supported in SLiM, I expect people will do models with even more traits, and they will want that to work with tree-sequence recording. So we clearly need a new design.

What we need is a proper place to put mutation metadata: one that can support binary data of arbitrary form, that can be written out at the end of a run (capturing current state) rather than at the moment a new derived state is created (capturing state that might be stale later on), and that can capture each mutation once rather than glomming things up with duplicate copies (as the derived state column does, for SLiM with stacking). And the question is where that table of metadata can live, in a .trees file, and how we write it out and read it in.

We could approach this as a general problem – maybe it makes sense for tskit to provide a general strategy for attaching metadata to mutations (without having to wedge it in as ASCII data in the derived state column), or maybe it makes sense for tskit to provide a general strategy for client code to add their own tables to a .trees file, and read those tables back in. Or we could approach this as a SLiM-specific problem – maybe this is really only SLiM's problem and nobody else's, in which case a general tskit-level solution is perhaps not needed, so then maybe we just figure out how to work at the kastore level on our own and write/read our own table behind tskit's back. Or maybe we even write out a separate file – a .slimtrees file – side by side with the .trees file, such that the pair of files are associated and the pair is needed to correctly read anything useful back in.

I don't know which avenue is better. To me, at least some kind of minimal sort of support in tskit for custom client tables seems like it would make sense, to make the .trees data structure extensible. That will probably prove useful to lots of people in lots of ways, going forward. But if that is a no-go for some reason, then we need to figure out how SLiM is going to hack this in without tskit support, I guess.

So, let's talk about this. I've raised this issue at points in the past, but it has seemed like a distant problem and has not gotten solved; now it has become rather urgent, since I'm doing the multitrait work in SLiM now and I need a solution. Seems like this is of interest to @petrelharp @jeromekelleher @benjeffery @hyanwong and probably others. I wasn't sure whether to make this issue on SLiM or on tskit; feel free to make a companion issue on the tskit side if it seems like that makes sense, but to avoid confusion let's either discuss here or there, not both. :->

[hyanwong]

Why doesn't SLiM use the actual metadata field for the mutations (or I guess the sites), rather than the derived state field? Is that because there is no easy C interface for metadata access? ISTR that @benjeffery has made it easy now to store binary ("struct") metadata in these fields and pull it all out (semi) efficiently for all tables (well, maybe not edges).

[bhaller]

Why doesn't SLiM use the actual metadata field for the mutations (or I guess the sites), rather than the derived state field? Is that because there is no easy C interface for metadata access? ISTR that @benjeffery has made it easy now to store binary ("struct") metadata in these fields and pull it all out (semi) efficiently.

Well, a site is different from a mutation, of course; even without the stacking of mutations that can happen in SLiM, there could be several different mutations segregating at a given site, each with different metadata. Now, I suppose we could keep metadata for all of the mutations at a particular site in the metadata entry for that site, and append new metadata to that site's metadata if a new mutation arose at that site. That might be better than using the derived state to keep all this metadata, since the site metadata is allowed to be binary rather than ASCII. It seems like it would still be a bit like pounding a round peg into a square hole; appending new metadata in the middle of the tskit metadata buffer would require quite a lot of copying of data and re-calculating offsets and so forth, so there would be tremendous overhead to maintaining it at runtime. Maybe we could keep all the metadata on the side in our own format, and then create the site metadata for tskit at save time. Anyhow, it's food for thought.

[hyanwong]

Well, you have the choice between saving in the mutation metadata field, or the site metadata field: perhaps a combination of both would work? I'm sure @benjeffery will have thoughts about the efficient way to write it (I don't!). I'm not sure if he is subscribed to messages in this repo though?

[bhaller]

Well, you have the choice between saving in the mutation metadata field, or the site metadata field: perhaps a combination of both would work?

The mutation table is not really well-suited to the task either, as described above. Fundamentally, that's because SLiM and tskit think about mutations differently, because SLiM allows stacking. An entry in the tskit mutation table therefore doesn't represent a SLiM mutation, but a whole SLiM stacked state; and that leads to the metadata for multiple mutations being recorded in a single derived state entry, and duplication of mutation metadata, and so forth. Neither the mutation table nor the site table is well-suited to SLiM's way of thinking about mutations, which is why we've had our hacky and unsatisfactory derived-state approach for so long.

I just want a place where SLiM can put its own mutation metadata. To me, the ideal would be if the derived state in SLiM simply recorded mutation IDs (for each mutation in the new stacked state), and then I had a place to put my own binary table with mutation IDs and their associated metadata. In particular, a way to add a new kastore table to a .trees file, used by the client code to store their own binary state without any involvement from tskit, seems like it would be ideal – and I doubt that SLiM will be the only client code that will want to do this, over time.

I'm sure @benjeffery will have thoughts about the efficient way to write it (I don't!). I'm not sure if he is subscribed to messages in this repo though?

I think tagging somebody on an issue flags it for them, regardless of where it is located? But I'll ping him on Slack just to make sure. :->

[hyanwong]

The tszip format allows arbitrary extra data to be stored in the file, and that's what is done in tsbrowse, ISTR. That might be an easy win, but it could be that you don't want to compress and uncompressed on the fly. Again, @jeromekellher and @benjeffery spent a long time figuring out some solutions here, so definitely worth checking to see if anything fits the bill.

[bhaller]

The tszip format allows arbitrary extra data to be stored in the file, and that's what is done in tsbrowse, ISTR. That might be an easy win, but it could be that you don't want to compress and uncompressed on the fly.

Interesting! I'm not sure whether that would be OK or not; I guess it depends on how fast the compression is? And whether there's any downside to the data being compressed, from the user's perspective?

Again, @jeromekellher and @benjeffery spent a long time figuring out some solutions here, so definitely worth checking to see if anything fits the bill.

Yep, for sure! That's what this chat is for! :-> Thanks @hyanwong, good to think about these various options!

[bhaller]

The tszip format allows arbitrary extra data to be stored in the file, and that's what is done in tsbrowse, ISTR. That might be an easy win, but it could be that you don't want to compress and uncompressed on the fly.

Interesting! I'm not sure whether that would be OK or not; I guess it depends on how fast the compression is? And whether there's any downside to the data being compressed, from the user's perspective?

I guess the other question, with respect to this, is what it would take for SLiM to write and read tszip files. Is that C/C++ code that I can throw into SLiM? Or is it Python-only, or a necessary external dependency, or some such? I have never looked into tszip at all. Homework for me, I guess, but right at this moment I need to go offline, as something else is happening now. :->

[benjeffery]

tszip isn't a valid route, there's no C API, and the extra tables will be ignored on decompress.

One option (off the top of my head) would be for the tskit C API to allow arbitary table data to go to the kastore, something like:

int tsk_table_collection_dump_with_extras(const tsk_table_collection_t *self, const char *filename, tsk_flags_t options, const tsk_kv_item_t *items, size_t num_items);

int tsk_table_collection_load_with_extras(tsk_table_collection_t *self, const char *filename, tsk_flags_t options, tsk_extras_reader_t reader, void *user_data);

If we wanted to extend to Python then it is something like:
ts = tskit.load(p, extras_processor=my_func); ts.dump(p2, extras=extras)

The standard tskit C API then transparently round-trips these extras. In the store I suggest paths like extensions/slim/mutations/blocks.

I'm not a massive fan of this as it invites all sorts of incompatible stuff to be built. If we can stick to standard mutation metadata it would be better - but I haven't fully understood why that doesn't work.

[petrelharp]

Some more background here: the basic issue is that in SLiM you can have complicated things like "mutations of type A and B can stack but do not stack with C" so that if an A happens on a B background then both are inherited, but then if a C occurs it removes the A and B so the offspring just gets the C mutation. Still, the allelic state of an individual would be inferrable from the history of mutations on the tree and the rules for stacking... except that the rules for stacking can change with time.

So, the cleaner thing would be, as @bhaller says, to have a new table that stores each unique new mutation (call this the "allele table??), since then each thing only occurs there once, and then in the mutation table what gets recorded in metadata is which rows in the allele table individuals inheriting that mutation (and not any subsequent ones) carry.

Currently we put a comma-separated list of SLiM mutation IDs as the derived state. If we wanted to avoid actually parsing this then we could stick the same information in metadata. The reason this is a good choice for derived state is that this string uniquely determines which mutations are carried, so if two are identical, then they are identical-by-state. We could make up some other system for this, I suppose.

This doesn't solve all the problems: as @bhaller says, information about mutations might change with time. For instance, what if all selection coefficients were randomly resampled every generation. However - as this example shows - we can't hope to come up with a one-size-fits-all solution that deals with this in a reasonable way. In designing the current system we decided (in many cases) stick to trying to only represent a 'snapshot' on SLiM objects; if they change with time then the user needs to store that some other way (for another example: what if an individual's spatial location changes with time).

[petrelharp]

AFAICT our options are:

1. keep with the current system, shove it all in mutation metadata
2. write out the per-mutation information to a separate table in the kastore somehow
3. write it separately (ick) or in top-level metadata (ascii, ick)

The downside of (1) is basically that the same information (about a given mutation) gets written down more than once (in the mutation metadata of any other mutation that appears on top of it); also, this is confusing for the end user. It's not clear to me how much this actually affects things, since at sensible mutation rates, multiple mutations at the same site are relatively rare. Even if 1% of sites have multiple hits, this is only going to increase memory usage for the mutation metadata by about 1%.

The downside of (2) is that it would mean having a whole new table. Also very confusing, and a lot of fuss to get working, and to describe in the docs, etcetera.

Looking again at what you say the problems are, @bhaller, the other problem with the current system is that "things might change with time". However, as in the previous message, I suspect we just can't reasonably deal properly with this. Are you thinking that if we had a new table (again, the 'allele table') that every time the attributes of a mutation changed, we'd write a new line in that table? So we'd need a time attribute in that table? That's fine in principle, but it'd be very easy for someone to blow up the table by changing things a lot. And we don't record time-changing things for any other objects (since it fundamentally breaks tskit's data model).

[petrelharp]

A different suggestion: what if mutation tables had a column that worked just like location does for individuals? Then SLiM could record the net effect of the various stacked mutations (if this is well-defined) on the various phenotypes, making things easier for the downstream user? I'm not sure this is worth it, but I'm brainstorming.

Also note that SLiM could start using more columns in location to record the phenotype values, which makes parsing the mutations less important.

[bhaller]

tszip isn't a valid route, there's no C API, and the extra tables will be ignored on decompress.

Aha, OK, good to know.

One option (off the top of my head) would be for the tskit C API to allow arbitary table data to go to the kastore, something like:

int tsk_table_collection_dump_with_extras(const tsk_table_collection_t *self, const char *filename, tsk_flags_t options, const tsk_kv_item_t *items, size_t num_items);

int tsk_table_collection_load_with_extras(tsk_table_collection_t *self, const char *filename, tsk_flags_t options, tsk_extras_reader_t reader, void *user_data);

If we wanted to extend to Python then it is something like: ts = tskit.load(p, extras_processor=my_func); ts.dump(p2, extras=extras)

The standard tskit C API then transparently round-trips these extras. In the store I suggest paths like extensions/slim/mutations/blocks.

Yes, I think this sounds like a reasonable solution. Of course details remain to be ironed out – tsk_kv_item_t, tsk_extras_reader_t – but I think I see what you're sketching here and it seems good. :->

I'm not a massive fan of this as it invites all sorts of incompatible stuff to be built.

I certainly see the objection; it opens the barn doors wide. But really the barn doors are already wide; people can put arbitrary stuff into a companion file, or into top-level metadata, etc. All that this proposal is really doing is providing a way of attaching such data in a form that is easier to work with – can be binary, doesn't need to be described by a JSON schema, doesn't need to be wedged into the derived state column or similar place, etc. – while still having it be packaged nicely inside the .trees file.

If we can stick to standard mutation metadata it would be better - but I haven't fully understood why that doesn't work.

Hopefully @petrelharp has helped above? Happy to zoom about it, or talk about it in a tskit meeting (but I really don't think it's of general interest to the group).