clarify derived state / metadata[derived states] distinction for mutations

[petrelharp]

Thanks to #665, now the SLiM IDs of mutations carried by a tskit mutation are recorded in two places: the mutation's derived_state, and in it's metdata["derived_states"]. Ideally, when reading things back in, we would ignore the former and only use the latter, but this is complicated because we're using the tsk_variant machinery, which returns alleles. So to use the metadata instead we'd need to either extend the tsk_variant behavior or write our own. This would be fine but also not top priority.

So currently, to read in the file we need the derived_state column to have the comma-separated list of mutations IDs in it. Are there any requirements on the metadata? Does the code check that these agree? We should record this in treerec/implementation.md and in the manual if it's not there already (I'm only partway through checking).

[bhaller]

The metadata is now documented as matching the derived state, but in binary rather than ASCII. Yes, this should be in implementation.md and should be checked by a python test; I can do the former, you should do the latter. :-> Sound good?

[bhaller]

OK, I just pushed a first pass at the implementation.md doc for this. It is already documented in the manual, I think; should be in the latest draft that I posted in stdpopsim. :->

[petrelharp]

From the docs:

I think this should say something like:

When SLiM reads in a tree sequence, it uses the derived state (not the metadata) to assign mutations to haplosomes.

and then perhaps,

It does not at present do anything with the list of derived states in mutation metadata (or even check for agreement between the list of mutation IDs in derived state and in metadata).

or the analogous statement if in fact it does do that consistency check?
Finally,

However, this behavior may change in the future, so both derived state and metadata should be present.

[petrelharp]

(This is kinda making me want to change the code so it uses metadata and not the derived state.)

[petrelharp]

The note in implementation.md in the recent commit is good (although I think a bit redundant with elsewhere in the file?) but it's missing the key note about what is the source of truth? (ie there's redundancy; so which is the source of truth and is agreement between them required?)

[bhaller]

I think this should say something like:

When SLiM reads in a tree sequence, it uses the derived state (not the metadata) to assign mutations to haplosomes.

and then perhaps,

It does not at present do anything with the list of derived states in mutation metadata (or even check for agreement between the list of mutation IDs in derived state and in metadata).

or the analogous statement if in fact it does do that consistency check? Finally,

However, this behavior may change in the future, so both derived state and metadata should be present.

It's the opposite, actually; so when you then write:

(This is kinda making me want to change the code so it uses metadata and not the derived state.)

That's what SLiM does. :-> See Species::DerivedStatesFromMetadata(). It ignores the ASCII derived state stuff (and doesn't check it for consistency, partly because that would be slow – eww, ASCII – and partly because it might've been converted to nucleotides or something – the user probably shouldn't do that, but I don't see a reason to actively block it). As far as SLiM is concerned, the binary derived state data is the truth.

(It's what I always really wanted the derived state info to be, but tskit doesn't allow binary! Actually, now that we have schemas, a nice idea would be to allow a derived-state schema to be specified in the tree sequence. SLiM could specify an "array of binary ints" schema, and then tskit would know how to print derived states from SLiM. Then maybe we could actually get rid of the ASCII derived states entirely! It'd make writing .trees files from SLiM a bit faster, and make the files a bit smaller, but more importantly it would just clean up all of this gunk about ASCII versus binary, derived state column versus metadata column, etc. Is this an idea worth considering?)

So, OK, I have added a new paragraph below the doc that you snapped above:

When SLiM reads in a tree sequence, it uses the binary derived-state information in the metadata column (not the derived state column) to assign mutations to haplosomes. At present SLiM ignores (and indeed, disposes of) the ASCII information in the derived state column, without checking for agreement between it and the binary derived-state information in the metadata column. However, this behavior may change in the future, so both derived state and metadata should be present, and should be in exact agreement.

Good?

[bhaller]

The note in implementation.md in the recent commit is good (although I think a bit redundant with elsewhere in the file?) but it's missing the key note about what is the source of truth? (ie there's redundancy; so which is the source of truth and is agreement between them required?)

Ah, I didn't realize that was what you were driving at. OK, I have just added this to implementation.md. See this commit.

[petrelharp]

That's what SLiM does. :-> See Species::DerivedStatesFromMetadata(). It ignores the ASCII derived state stuff (and doesn't check it for consistency, partly because that would be slow – eww, ASCII – and partly because it might've been converted to nucleotides or something – the user probably shouldn't do that, but I don't see a reason to actively block it). As far as SLiM is concerned, the binary derived state data is the truth.

Well gee, I thought SLiM used metadata and ignored ASCII, glad I wasn't totally wrong.

But let's see: tsk_variant_t appears in the following places in species.cpp:

• in CrosscheckTreeSeqIntegrity
• in TallyMutationReferencesWithTreeSequence
• in AddMutationsFromTreeSequenceToHaplosomes

And, the tsk "variant" stuff grabs site ancestral states and mutation derived states to do its thing.

I noticed this because if I do pyslim.convert_alleles to put nucleotides in the derived state column (and into site ancestral states) and then try to load the tree sequence into SLiM, I hit the associated error:

E           ERROR (Species::__TallyMutationReferencesWithTreeSequence): (internal error) variant allele had length that was not a multiple of sizeof(slim_mutationid_t).

Former us realized this, I think: see the second half of this comment

So, I don't think I agree with your additions above?

But, I may well be missing something?

[petrelharp]

Then maybe we could actually get rid of the ASCII derived states entirely! It'd make writing .trees files from SLiM a bit faster, and make the files a bit smaller, but more importantly it would just clean up all of this gunk about ASCII versus binary, derived state column versus metadata column, etc. Is this an idea worth considering?

Very tempting, but the current scheme has the very nice attribute that individuals have the same allelic state (which is determined by the derived_state if and only if they have the same collection of SLiM mutations. If we changed this, we'd get different results for, for instance, popgen statistics calculated from the tree sequence. Other schemes besides "comma-separated list" would of course be possible but why overcomplicate things.

[bhaller]

That's what SLiM does. :-> See Species::DerivedStatesFromMetadata(). It ignores the ASCII derived state stuff (and doesn't check it for consistency, partly because that would be slow – eww, ASCII – and partly because it might've been converted to nucleotides or something – the user probably shouldn't do that, but I don't see a reason to actively block it). As far as SLiM is concerned, the binary derived state data is the truth.
...

But let's see: tsk_variant_t appears in the following places in species.cpp:
...

SLiM has always kept binary derived state info in the derived state column and converted it to ASCII at save time, since we first developed tree-sequence recording in SLiM. The variant stuff doesn't mind, I guess; it just uses the binary bits in the derived state column and doesn't care what they represent. The only place where tskit cares that the derived state info is ASCII, as I understand it, is when printing derived states, and that doesn't arise while SLiM is running. So keeping binary in the derived state column is safe.

I noticed this because if I do pyslim.convert_alleles to put nucleotides in the derived state column (and into site ancestral states) and then try to load the tree sequence into SLiM, I hit the associated error:

E           ERROR (Species::__TallyMutationReferencesWithTreeSequence): (internal error) variant allele had length that was not a multiple of sizeof(slim_mutationid_t).

Former us realized this, I think: see the second half of this comment

So, I don't think I agree with your additions above?

But, I may well be missing something?

I'll look into that. But I think it'll just be a bug of some sort – fallout from the recent changes to the way the derived state information is shuffled around – not anything significant, or relevant to this discussion.

[bhaller]

Then maybe we could actually get rid of the ASCII derived states entirely! It'd make writing .trees files from SLiM a bit faster, and make the files a bit smaller, but more importantly it would just clean up all of this gunk about ASCII versus binary, derived state column versus metadata column, etc. Is this an idea worth considering?

Very tempting, but the current scheme has the very nice attribute that individuals have the same allelic state (which is determined by the derived_state if and only if they have the same collection of SLiM mutations. If we changed this, we'd get different results for, for instance, popgen statistics calculated from the tree sequence. Other schemes besides "comma-separated list" would of course be possible but why overcomplicate things.

No, I'm just proposing that the derived states could be binary instead of ASCII, if there was a schema attached to the derived state column. Otherwise the actual information content in the derived state column would be unchanged. Just binary instead of ASCII. But printable, thanks to the schema.

[bhaller]

I noticed this because if I do pyslim.convert_alleles to put nucleotides in the derived state column (and into site ancestral states) and then try to load the tree sequence into SLiM, I hit the associated error:
...

Can you post the tree sequence and the SLiM model you're loading it into, so I can reproduce? Thanks!

[petrelharp]

Okay, gotcha. I was only talking about the input/output tree sequences, not the internal stuff, and the reason was because of this issue I ran into. Here's the script:

initialize() {
	defineConstant("L", 1e2);
	initializeSLiMOptions(keepPedigrees=T, nucleotideBased=T);
   initializeTreeSeq(timeUnit="generations");
	initializeAncestralNucleotides(randomNucleotides(L));
	initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
	initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(4e-2));
	initializeGenomicElement(g1, 0, L-1);
	initializeRecombinationRate(1e-2);
}


1 first() { // reference sequences
	// sim.readFromPopulationFile("pre_cts.txt"); // succeeds
	sim.readFromPopulationFile("cts.txt");  // fails
}


20 late() {
    catn("Done.");
}

and here are the tree sequences:
this one works: pre-cts.txt
and this one fails: cts.txt
with

ERROR (Species::__TallyMutationReferencesWithTreeSequence): (internal error) variant allele had length that was not a multiple of sizeof(slim_mutationid_t).

[bhaller]

Okay, gotcha. I was only talking about the input/output tree sequences, not the internal stuff, and the reason was because of this issue I ran into. Here's the script:
...

ERROR (Species::__TallyMutationReferencesWithTreeSequence): (internal error) variant allele had length that was not a multiple of sizeof(slim_mutationid_t).

OK, hmm. Reproduces fine on my side. The mutation table is getting cleaned up correctly; after DerivedStatesFromMetadata() it correctly has the binary derived state info in the derived state column, and the metadata column is empty. That setup is how things have always been for SLiM while it is running. The "variant allele had length that was not a multiple of sizeof(slim_mutationid_t)" error appears to be coming from somewhere else – maybe tsk_variant_decode() is picking up an ASCII allelic state from the site table, not the mutation table? Does pyslim.convert_alleles() do something there? I imagine it does, and that that is what is causing the problem...?

In the past we haven't supported loading a .trees file back in after convert_alleles(), right? So this is a new thing we're trying to do. Is it worthwhile to do it, or should we just say – as we have before – "you can't do that"? If we want to do it, then I guess we need to decide what to do with the site table. Would SLiM want to just clear out the derived state column there, making all the ancestral states zero-length? I guess that is the state that SLiM expects, and it's certainly not going to understand/use nucleotides that are in there.

Unless we've got a compelling reason to support this new workflow, I think we can just punt on it. What do you think? The only thing that might be nice would be to improve the error message, I guess; if the unexpected variant length is == 1, the error could add "(If you called pyslim.convert_alleles() on the tree sequence before loading it in to SLiM, note that that is not supported because it leaves the tree sequence in a state that SLiM cannot parse.)"

[petrelharp]

Ah right: it could be the site table - convert_alleles puts a nucleotide in the ancestral state there, and that'll show up in the tsk_variant stuff.

I think this is probably a very easy fix - just remove the ancestral states of the sites on loading the tree sequence in. I'll check it out.