Ancestral states when sample data is missing

[hyanwong]

Something I just realised, after implementing a unit test for tsinfer. If we have samples that are missing a left hand or right hand chunk of information, then my recent push (tskit-dev/tsinfer#169) allows us to build a decent tree sequence, and a samples that has a particular region missing will simply have its node number unconnected to any other node in that region of the tree sequence (it will just be an unconnected point when plotted). But when we generate a haplotype from this tree sequence, there will be zeroes not tskit.MISSING_DATA for the missing regions. That’s because the ancestral state is zero in these regions, and since they aren’t linked by an edge to anywhere, they just get the ancestral state.

There are a few reasonable ways around this. One way (1) is to say that any completely unconnected nodes should always get a missing data flag. Another (2) is to say that for these sites, the ancestral state is tskit.MISSING_DATA, then place a mutation to ‘0’ above the root of the main tree. Yet another (3), which is more generic but has a wider impact is to set the ancestral state to tskit.MISSING_DATA for all sites, and always require a mutation to 0 to be placed above the main tree root.

I see this as a general problem that if we have a position in the tree sequence where there are multiple roots, there might be different ancestral states for each root.

[jeromekelleher]

I see this as a general problem that if we have a position in the tree sequence where there are multiple roots, there might be different ancestral states for each root.

Well, there would still be an ancestral state, it's just that some of the roots would have mutations over them. The concept of a state that is ancestral to all of the samples is still a sound one, it's just that we're missing some of the topology.

This is only an issue when we're generating genotypes for sample nodes that have incomplete data, and we could add this behaviour with an appropriate flag to variants.

[hyanwong]

I guess I'm unsure if by "ancestral state" we mean the state ancestral to all the samples (as you say), or of only the samples for which we have data. The first is essentially unknowable. The second is not. But if you want the definition to be the first one, then I guess the most principled thing would indeed be to set the "ancestral state" to "UNKNOWN", and then simply place a mutation over the root of the known samples somehow? Or have you another idea?

we could add this behaviour with an appropriate flag to variants.

Do you mean on input, or on output? I'm getting the matrix out by doing ts.genotype_matrix(), and it would be nice if this behaved properly, so that I can round-trip the data, passing missing data into a SampleData instance, and getting the missing states back out from the resulting inferred ts.

[hyanwong]

@jeromekelleher - I guess the question here is do we:

a. want to encode missing data in a tree sequence by using the existing format (e.g. setting ancestral_state to tskit.MISSING_DATA and placing a mutation to 0 above the "main" root at that site), OR

b. do we want to alter tskit in some way, e.g. for isolated nodes with no edges connecting them to a tree at a given location, output tskit.MISSING_DATA. or perhaps by passing something to the ts.variants() method.

Option a. can be done by altering tsinfer. The options under b. require adjustments to tskit, I presume.

[jeromekelleher]

(a) is an elegant solution @hyanwong, I didn't think of this! I don't think we can change the current encoding of sites and mutations without breaking code, so we should be very careful about doing it. I guess it's reasonable thing to do that if you are in a situation where missing data is an issue that you encode this as the ancestral state of the site. For simulations, we always know the ancestral state so we can just encode it directly.

Do we need to make any changes in tskit to get this working? I'm not sure we'll allow it at the moment as the ancestral state must be a string. We'd have to think carefully how to enable it, and make sure we're not breaking code.

The simplest thing we can do to allow us round-trip genotypes is to add a flag to the variants method, which is what we should do first I think.

[hyanwong]

I'm not sure we'll allow it at the moment as the ancestral state must be a string.

That's a good point. I presume that for e.g. indels we might want to have empty strings as one of the alleles, so we can't use an empty string to mark missing data. I also assume that allele info is stored in some sort of zarr array, so we can't use None or some other python indicator. A semi-nasty hack might be to define tskit.MISSING_ALLELE_STATE = chr(127), or some other ASCII value that won't ever be used in an allele string and which could reasonably be thought of as "missing", such as ESC (chr(127)), EOT (chr(127)) or CAN (chr(24)).

The simplest thing we can do to allow us round-trip genotypes is to add a flag to the variants method, which is what we should do first I think.

So how about adding a parameter for a dictionary that maps strings to "alternative" genotype numbers? E.g. variants(alternative_mapping = {chr(127):-1}), meaning that any genotype value mapping to the string chr(127) should be given the genotype value -1, rather than the normal index into the genotype array.

There is a decision to be made here, about whether we keep all the other values the same in the genotype array. If the ancestral state was an allele of chr(127), then we would also have two other alleles (e.g. A and T) with genotype values 1 and 2, so if we keep all other values the same, we would end up with variants with no 0 values, e.g.

(-1, -1, 1, 2, 2, 2, 1, -1)

rather than the original sample data file which would be

(-1, -1, 0, 1, 1, 1, 0, -1)

Alternatively we could re-index the alleles, removing any marked as missing, but that would mean altering the alleles list, which seems horrible to me.

Do any other solutions occur to you?

[jeromekelleher]

I don't think we need a string corresponding to MISSING_DATA if we add a flag to variants - we just return -1 and either add None/NULL to the end of the alleles list or leave it as it is. I'll code up a proposal when I get a chance, should be straightforward to do and it'll be really useful for tsinfer.

[hyanwong]

Thanks for offering to code up, as I don't quite understand your proposal! I can't see how we would encode missing data in the tree sequence itself, without reference back to the original sample data file, so I can't see how this would properly round-trip, but I think I'm simply being dense. I'll wait to see your proposal.