Need a flag to decide what to do about sites that are not shared between multiple VCZ input files

[hyanwong]

It will be common to have (say) 2 separate VCZ files, A and B with overlapping but not identical sites. For sites that are in A.vcf but not in B.vcf, we need to decide what to do. The options are basically to treat the allelic state for the samples in B.vcf as missing, to treat them as REF, or to treat them as the ancestral state.

At the moment, the design.md file has a draft for phase 2 that states "Genotypes at sites absent from a source are treated as missing for that source's samples." Although missing seems reasonable (and is probably what we did for the unified genealogy), I suspect this will lead to an incorrect frequency value, because the most likely interpretation is that the site is entirely REF in B.vcf. This is a different sort of missingness to case where we actually don't know the genotype for specific samples in sites defined within a VCF.

So I suggest that we update the design and mark these sites as REF by default. We could either have a parameter that allows this to be changed so that they are treated as missing, or we could simply require the use to define an all-missing site if they want this behaviour.

As an aside, we should also check that the sites that are shared between the VCFs define the same REF, and raise an error if not. They need not have the same ALT alleles, however.

[hyanwong]

Alternatively, a more sophisticated approach would be to mark these as missing in the genotype matrix, but within tsinfer, calculate the frequency (used to time order the ancestral haplotypes) by assuming that the missing data values are the REF.

[hyanwong]

I think that we did this wrongly in the unified genealogy. For reference, the missingness code that @awohns used is at https://github.com/awohns/unified_genealogy_paper/blob/f71155e2b8b2cc096607803e563de48d215b4c29/all-data/Makefile#L671, and a bit of sleuthing shows that links to

tsinfer/tsinfer/formats.py

Line 1951 in de57136

def merge(self, other, **kwargs):

. From this, it seems to me that sites in one dataset but not the other are treated as missing, which as argued above,, is probably not correct.

[hyanwong]

The reduction ad absurdum for this is where each of N diploid samples is in a separate VCF, and we merge all N VCFs together. In this case, we definitely don't want to treat a heterozygous site present in only one VCF as present at 50% frequency in the population as a whole. This would be a good unit test case.

[Duncan-JR]

I discussed with @hyanwong and we agreed on what we think should be the ideal behaviour.
If $D$ is the number of derived alleles of a variant across all sources, $N_{present,nonmissing}$ is the total number of non-missing genotypes in sources containing the variant, and $N_{absent}$ is the number of genotypes in sources where the variant is absent, the focal frequency should be calculated as

$$ f = \frac{D}{N_{present,nonmissing} + N_{absent}} $$

For ancestor building, we agree that sources where a variant is missing should have missing genotypes for that variant, which is the current behaviour. Since our discussion, I discovered that the current build simply calculates frequencies as total derived AC / total number of genotypes across all sites. This means that frequencies are calculated correctly with N diploids in separate VCFs, but actual missing genotypes are always included in the calculation as well (unlike in tsinfer 0.5).

There is a caveat to consider: since we compute the derived allele frequency, if a site is absent from a source, we have to check whether the REF allele is ancestral or not. If it isn't, I think we should count the total number of per-source genotypes in the numerator as well as denominator (since we assume all are REF for the purposes of allele frequencies). So, when the multi-source view is constructed for ancestor generation, I think we need to compute:

$$ f = \frac{D + N_{absent,ref is derived}}{N_{present,nonmissing} + N_{absent}} $$

@hyanwong does this look good to you? Implementing it shouldn't be too tricky. I suspect the most time-consuming part will be setting test cases.