Add an "ancestral_state_index" parameter to SampleData.add_sites()?

[hyanwong]

At the moment we require people to order their allele lists such that the first allele is the ancestral state. This can be a little fiddly, and often the ancestral state isn't known anyway, so it would be good to have a way for users to indicate this, while keeping the site in the final TS. Currently the way to incorporate a site with an unknown ancestral state is to pick one at random as the first one to put, and include the site in the exclude_positions parameter of generate_ancestors(), which seems pretty obscure to me.

I suggest that a non-breaking change that could allow this would be to have an extra parameter "ancestral_state_index", which defaults to 0, meaning that the zeroth allele in the alleles parameter is treated as the ancestral state. This would be similar to the sort of format we expect to see when we switch to using sgkit as the input format. We could allow tskit.NULL to be passed to mean that we don't know the ancestral state, and it should be placed using parsimony.

This seems pretty intuitive to me, but there are some decisions to be made here. When we add a site with (say) ancestral_state_index=1, do we record the genotypes and allele list before sticking them into the sample data file (which seems a bit hacky, and is not likely to be portable to sgkit). Or do we switch the genotypes when we build the ancestors, in which case we will need to swap the genotypes when we match against the ancestors using match_samples or match_ancestors. If the latter, I suspect that it might be helpful for sanity checking purposes to store the meaning of the genotype values (i.e. the reordered alleles list) in the AncestorData instance.

[jeromekelleher]

Nice idea, I like it. I think we should separate out the idea of saving people from fiddling with the allele ordering from specifying whether the ancestral state is known or not though.

Or, we could make it simpler yet, by giving the actual allele value of the ancestral state, so that people don't have to figure out what index its at at all?

Ultimately we're probably going to want to have an ancestral_state_confidence or something, and we do something slightly different for low, medium and high-confidence ancestral state calls (which is why I don't think we should mush the functionality into this proposal).

[hyanwong]

Or, we could make it simpler yet, by giving the actual allele value of the ancestral state, so that people don't have to figure out what index its at at all?

Yes, that works. None (default) would mean take the first in the allele list. But although I agree with your general point about separating out functionality, I do think we should consider allowing tskit.NULL here, as it's weird to have to specify an ancestral state if we really don't know it at all.

Ultimately we're probably going to want to have an ancestral_state_confidence or something, and we do something slightly different for low, medium and high-confidence ancestral state calls (which is why I don't think we should mush the functionality into this proposal).

I agree with that too, and it's worth bearing this in mind. I think that would probably be added to both the SampleData and the AncestorsData formats.

[jeromekelleher]

I do think we should consider allowing tskit.NULL here, as it's weird to have to specify an ancestral state if we really don't know it at all.

We can't use tskit.NULL as that's -1, which people could quite reasonably use to mean the last element of the array.

I agree with that too, and it's worth bearing this in mind. I think that would probably be added to both the SampleData and the AncestorsData formats.

I think we'd end up having to change these anyway to do what you're talking about, so lets do the simple thing first.

[hyanwong]

OK, so do we store the ancestral state (and presumably the confidence, when we implement it) in both the SampleData and AncestorsData files? Or do we simply reorder the alleles and recode the genotypes in the SampleData file?

If we go down this route, I'm slightly inclined to do the confidence thing first, TBH, as we have a reasonable workaround for the order thing, but not for adding sites where the ancestral state is unknown. I guess we would store the confidence as a float (0..1)

[hyanwong]

I do think we should consider allowing tskit.NULL here, as it's weird to have to specify an ancestral state if we really don't know it at all.

We can't use tskit.NULL as that's -1, which people could quite reasonably use to mean the last element of the array.

Hmm, I just thought: tskit.NULL would be fine if we expected a value ("A", "T", "G", "C"), rather than an index, right?

[hyanwong]

Summarising the discussion above, there seems like a reasonably easy API for this: we add an ancestral_state parameter i.e. sd.add_site(..., alleles=..., ancestral_state=None) where ancestral_state is a string that matches one of the alleles, or None to mean the zeroth allele, or tskit.NULL to mean we don't know. We leave the confidence stuff for a future PR. We add the ancestral_state as a new column in the sample data file (or if we only want to add an int column, we store the index instead).

The gnarly bit is what then to do with it. If there's a simple way to integrate this information into the current ancestor generation / matching algorithm, then I think we could reasonably add it to https://github.com/tskit-dev/tsinfer/milestone/3. If it's not simple, I suggest punting it off to a later milestone.

[jeromekelleher]

After discussion, decision is to add the new column to the SampleData file, and to bump the minor version number. We can assume that files missing this column have an ancestral_state of all zeros.

[hyanwong]

I think in tsinfer we require the derived state to be 1, right? So I suggest that we request both an ancestral state and a derived state from the user (but default to using the zeroth and first allele in the list), then store both a column of ancestral states and a column of derived states in the sample data file. That way, as long as the other alleles are left in an unchanged order, we can output the same variant encoding as we passed in. For the same reason as discussed, this should only require a minor version number bump. Here's an example:

sd.add_site(position=0, genotypes=[2, 3, 2], alleles=["A", "T", "G", "C"], ancestral_state="C", derived_state="G")
sd.sites_genotypes[0] == [1, 0, 1]  # suitable for inference: 0 is the ancestral state
sd.sites_alleles[0] == ("C", "G", "A", "T")  # NB: "A" and "T" are in the same order as passed in
sd.sites_original_ancestral_allele_index[0] == 3  # defaults to 0 if `ancestral_state` not specified in add_site
sd.sites_original_derived_allele_index[0] == 2  # defaults to 1 if `ancestral_state` not specified in add_site

Then we can make sure that the variants() operator by default returns the same order as passed in, even though under the hood it is storing a different array:

def variant(self, site_id):
    # by default, return the genotype encoding that we passed in. We could also add a
    # param to return the "internal" encoding, e.g. for using in tsinfer, 
    # or another encoding via a passed-in allele list, like ts.variants()
    genotypes = sd.sites_genotypes[site_id]
    alleles = sd.sites_alleles[site_id]
    allele_index = [
        sd.sites_original_ancestral_allele_index[site_id],
        sd.sites_original_derived_allele_index[site_id],
    ]
    allele_map = np.array(allele_index +[i for i in range(len(alleles) if i not in allele_index])
    # allele_map = [3, 2, 0, 1] in the example above
    reverse_map = {v: k for k, v in enumerate(allele_map)}
    genotypes = allele_map[genotypes]
    alleles = tuple(alleles[reverse_map[i]] for i in range(len(alleles)))
    return self.site(site_id), alleles, genotypes

# sd.variant(0) == 0, [2, 3, 2], ("A", "T", "G", "C")

Shall I work this up into a PR at some point? I think it then is completely backwards compatible.

[jeromekelleher]

I don't see why we need the derived_state - the matching engine can handle arbitrary alleles, and the ancestor generation only works for biallelic sites?

I don't think we need to preserve the input variant encoding - there will be plenty of cases where the ancestral state is the ALT allele, e.g.

[hyanwong]

Does the matching engine work with sites arrays like [0, 3, 3, 3, 0, 0] then? I thought it required biallelic sites to encode the derived state as a 1?

[hyanwong]

I don't think we need to preserve the input variant encoding

Not preserving the original encoding was causing complete havoc when Shing was trying e.g. to merge sample data files, mess with tree sequences vs SD files, etc. I thought this might cause a lot less pain for end-users.

[jeromekelleher]

Does the matching engine work with sites arrays like [0, 3, 3, 3, 0, 0] then? I thought it required biallelic sites to encode the derived state as a 1?

I think that's fine, it just treats them as numbers to match. There might be some weirdness with how it decides what the number of alleles is though, come to think of it.

[hyanwong]

Does the matching engine work with sites arrays like [0, 3, 3, 3, 0, 0] then? I thought it required biallelic sites to encode the derived state as a 1?

I think that's fine, it just treats them as numbers to match. There might be some weirdness with how it decides what the number of alleles is though, come to think of it.

I think the ancestors generator assumes 0/1 though?

[jeromekelleher]

I think the ancestors generator assumes 0/1 though?

I don't know.