Tskit: Genotype encoding scheme

[jeromekelleher]

The implementation of the parsimony method over in #125 has raised some troublesome issues.

1. We need to be able to represent missing data
2. If we want to provide a conventional cost matrix for e.g., Sankoff parsimony, then genotype encoding 0, 1, 2... needs to have a fixed interpretation

Currently, the genotypes output by the variants method are encoded in a site-specific way, and the alleles array is needed to decode the actual string representation. A genotype of 0 at one site might mean "A" and at another is might mean "ATTAAC". This is fine and works well for producing data because (a) it gives us a lot of flexibility in terms of representing stuff like short indels and (b) always have the ancestral state encoded as 0, with other alleles numbered as you go down the tree is a useful property.

For inputting the data, I suggest we require the following to be provided:

• An input array of genotypes (np.uint8)
• A list of alleles
• For Sankoff parsimony, a len(alleles) x len(alleles) cost matrix (we can provide some higher level allele-> allele mapping, which would take the pain out of this? I.e., cost={"A": {"A": 0, "C": 0.25,...}...}

As I see it, there are then two options for specifying missing data: either -1 or len(alleles). Encoding as -1 seems nicer, but genotypes are currently encoded as uint8, so it would really be 255. Using 255 as the missing data value seems quite nasty to me, as we may need to add support for 16 bit genotypes at some point (already available in the C code), and suddenly 255 wouldn't mean missing data any more. This seems quite ugly, and may lead to tricky bugs. So, if we really do want to use -1 to represent we should change the type of genotypes to int8, so that we can represent -1 properly. This would mean that we only have space for 128 alleles, and may mean needing to have 16 bit genotypes sooner. This would be a bit of hassle to implement, and may break some people's code.

The other option is to use len(alleles) as the missing data value. The downside here is that you can't scan the genotypes for missing data, without knowing the alleles.

Any thoughts?

[petrelharp]

Clarification: when you say "for input" this is currently only for the parsimony algorithm, but would be for other things in the future?

[petrelharp]

Another option is to have "missing" be a special allele or just a value in the list of alleles. This isn't ridiculous because sometimes "missing" is an allele itself (e.g., genotype at a SNP for an individual with a large deletion there).

[jeromekelleher]

Clarification: when you say "for input" this is currently only for the parsimony algorithm, but would be for other things in the future?

Yes, currently for just the parsimony algorithm, but I'll be using it for the Li and Stephens haplotype matching code as well. I'd imagine there might be other situations when we might want to take haplotype data as input?

Another option is to have "missing" be a special allele or just a value in the list of alleles. This isn't ridiculous because sometimes "missing" is an allele itself (e.g., genotype at a SNP for an individual with a large deletion there).

We do need to know which allele it is though, as missing is treated specially. For example, in parsimony, it's treated as "can be any allele". This is useful where we want to impute data using the tree: if we have a dodgy site we're not sure of, we mark it as missing and its value gets imputed in a very natural and neat way.

[hyanwong]

I prefer having tskit.NULL as the missing data value. In fact, I think there are several sorts of missingness, e.g. "we know there is no value for the sequence at this site" (e.g. a SNP in an indel) vs "we have not sequenced this bit of the genome yet". I think there may be some others too. Allowing negative numbers means we have a few possibilities for different encodings for all these missing types.

I don't think restricting to 128 values is a problem. If users start hitting this, then they are probably in a weird regime where 256 would be not enough either.

When we do Li & Stevens matching in tsinfer, I think the ancestral haplotypes produced have -1 for missing regions, don't they, so this would also be consistent between tsinfer and tskit. In fact, I can imagine examples where the tsinfer and tskit notation needs to be the same (e.g. we want to do inference including ancient sequences, which have missing regions, and we want to pretend that these samples could be ancestors too).

How much would it break to make genotypes np.int8 vs np.uint8?

[hyanwong]

I think the ancestral haplotypes produced have -1 for missing regions, don't they,

I was wrong here. In tsinfer we currently use constants.UNKNOWN_ALLELE = 255, so my argument about tskit vs tsinfer consistency doesn't apply.

I still think that both tskit and tsinfer should have an explicit set of "missing" values for genotypes, and using negative numbers is probably the least confusing way to do this.

[petrelharp]

Seems like we just need a tskit.MISSING_ALLELE; and I think that setting this as -1 and changing to int8 seems fine.

[jeromekelleher]

OK, that's decided then: -1 is the value and we need to change the types for genotypes to ints. Shouldn't cause much problems I'd imagine.

[hyanwong]

Right. It might be useful to address https://github.com/tskit-dev/tskit/issues/192 first, which should help enforce consistent behaviour. I'll work it into a PR.