Mismatch between haplotype output and genotype output

[XingerTang]

AlphaImpute2 version: v0.0.2

With an input pedigree as follows:

1 0 0
2 0 0
3 1 2
4 1 3

and an input genotype as follows:

1 9 2 0 0 0
2 0 2 0 0 2
3 1 9 9 9 9
4 1 2 0 0 0

By running the following command

AlphaImpute2 -genotypes tests/functional_tests/test_2/genotypes.txt -pedigree tests/functional_tests/test_2/pedigree.txt -ped_only -phase_output -seg_output -out tests/functional_tests/outputs/test_2

We get the following genotype output:

1 1 2 0 0 0
2 0 2 0 0 2
3 1 2 0 0 1
4 1 2 0 0 0

and the following haplotype output:

1 1 1 0 0 0
1 1 1 0 0 0
2 0 1 0 0 1
2 0 1 0 0 1
3 1 1 0 0 0
3 0 1 0 0 1
4 1 1 0 0 0
4 1 1 0 0 0

We can see that the first loci of both the individual 1 and the individual 2 have their haplotypes and genotypes unmatched.

[gregorgorjanc]

Good catch @XingerTang!

Let's see what we could expect. For simplicity I have merged pedigree and genotypes

1 0 0 | 9 2 0 0 0
2 0 0 | 0 2 0 0 2
3 1 2 | 1 9 9 9 9
4 1 3 | 1 2 0 0 0

Phasing the homozygotes gives

1 0 0 | 9 1 0 0 0
1 0 0 | 9 1 0 0 0
2 0 0 | 0 1 0 0 1
2 0 0 | 0 1 0 0 1
3 1 2 | 0/1 9 9 9 9
3 1 2 | 0/1 9 9 9 9
4 1 3 | 0/1 1 0 0 0
4 1 3 | 0/1 1 0 0 0

Ignoring recombination one outcome for 3 could be

3 1 2 | 9 1 0 0 0 from 1.1 or 1.2
3 1 2 | 0 1 0 0 1 from 2.1 or 2.2

But since first genotype of 3 is 1, this suggests

3 1 2 | 1 1 0 0 0 from 1.1 or 1.2
3 1 2 | 0 1 0 0 1 from 2.1 or 2.2

and then propagating this upwards

1 0 0 | 1 1 0 0 0
1 0 0 | 1 1 0 0 0
3 1 2 | 1 1 0 0 0 from 1.1 or 1.2
3 1 2 | 0 1 0 0 1 from 2.1 or 2.2
4 1 3 | 0/1 1 0 0 0
4 1 3 | 0/1 1 0 0 0

Since both 1 and 3 are parents of 4 we could have

4 1 3 | 1 1 0 0 0 from 1.1 or 1.2
4 1 3 | 1 1 0 0 0 from 3.1

Taken together

1 0 0 | 1 1 0 0 0
1 0 0 | 1 1 0 0 0
2 0 0 | 0 1 0 0 1
2 0 0 | 0 1 0 0 1
3 1 2 | 1 1 0 0 0 from 1.1 or 1.2
3 1 2 | 0 1 0 0 1 from 2.1 or 2.2
4 1 3 | 1 1 0 0 0 from 1.1 or 1.2
4 1 3 | 1 1 0 0 0 from 3.1

Which matches the output that you see for haplotypes, but indeed the genotypes don't match these - 1st genotype of individuals 1 and 4 seems off! Summing up the above I get

1 0 0 | 2 2 0 0 0
2 0 0 | 0 2 0 0 2
3 1 2 | 1 2 0 0 1
4 1 3 | 2 2 0 0 0

While input genotypes are

1 0 0 | 9 2 0 0 0
2 0 0 | 0 2 0 0 2
3 1 2 | 1 9 9 9 9
4 1 3 | 1 2 0 0 0 first one clashes with the result (a genotype error?) 

I don't take multiple probabilities into account in the above. Is it possible that there are other possibilities and probabilities were such that genotype probability rounding caused an issue or there is just a plain bug in here?

[XingerTang]

@gregorgorjanc This is a result of AlphaImpute2's preference in handling cases that are heterozygous, where both haplotypes are likely to be 1. During the pedigree imputation, letting both haplotypes be the called state actually improves the overall accuracy of the outputs for large inputs, but may lead to this mismatching situation.
My fix is just a simple check of all the heterozygous genotypes after the computation to choose one of the haplotype to be 1 and the other one to be 0, to keep the outputs consistent.

[XingerTang]

The new commit 3c885be uses genotype probabilities to choose the most likely haplotypes.
To enable this, genotype probabilities for individuals with no offspring are added, which were previously missing.
The genotype results of the ped-only mode and the pop-only mode are not affected by the correction of the haplotypes.
While for the combined mode, these changes in the genotype probabilities negatively impacted the genotype accuracies for the early generations and positively impacted the genotype accuracies for the later generations. This might be because the inclusion of the genotype probabilities of the individuals without progenies is only going to have a strong correlation with the later individuals, but not the ancestors. Since the combined mode uses the outputs from the pedigree imputation to help the population imputation, the inclusion of these haplotypes may reduce its accuracy in the population imputation for the early generations.

Nevertheless, the differences are small and can be negligible.

Before:

[1] "outputs/ai2.genotypes"
     generation yield yieldMissing       acc
[1,]          0     1            1 0.9613021
[2,]          1     1            1 0.9944063
[3,]          2     1            1 0.9982353
[4,]          3     1            1 0.9974643
[5,]          4     1            1 0.9980814
[1] "outputs/pop_only.genotypes"
     generation yield yieldMissing       acc
[1,]          0     1            1 0.9608854
[2,]          1     1            1 0.9916563
[3,]          2     1            1 0.9937132
[4,]          3     1            1 0.9914643
[5,]          4     1            1 0.9915116
[1] "outputs/ped_only.genotypes"
     generation   yield yieldMissing       acc
[1,]          0 0.64045   0.06820313 1.0000000
[2,]          1 0.83213   0.47540625 1.0000000
[3,]          2 0.99614   0.98584559 0.9998881
[4,]          3 0.99654   0.98775000 0.9996023
[5,]          4 0.99437   0.98363372 0.9999704

After:

[1] "outputs/ai2.genotypes"
     generation yield yieldMissing       acc
[1,]          0     1            1 0.9594271
[2,]          1     1            1 0.9919062
[3,]          2     1            1 0.9982353
[4,]          3     1            1 0.9977857
[5,]          4     1            1 0.9983140
[1] "outputs/pop_only.genotypes"
     generation yield yieldMissing       acc
[1,]          0     1            1 0.9608854
[2,]          1     1            1 0.9916563
[3,]          2     1            1 0.9937132
[4,]          3     1            1 0.9914643
[5,]          4     1            1 0.9915116
[1] "outputs/ped_only.genotypes"
     generation   yield yieldMissing       acc
[1,]          0 0.64045   0.06820313 1.0000000
[2,]          1 0.54200   0.47540625 1.0000000
[3,]          2 0.99567   0.98654412 0.9998882
[4,]          3 0.99652   0.98810714 0.9996024
[5,]          4 0.99454   0.98485465 0.9999705

[gregorgorjanc]

@XingerTang the difference is indeed small. How does it impact the accuracy for haplotypes?

[XingerTang]

@gregorgorjanc I'm still using the AlphaImpute2 example data to test the accuracy. It only has the true genotype.
I think it should be better since fewer unexpected 1s are in the haplotype outputs.

[gregorgorjanc]

@XingerTang I guess this is now a good push to implement the more rigorous testing you implemented in AlphaPeel;) I do agree with your assessment and the code changes though! But having the test cases is kinda important, as you can see.

[XingerTang]

@gregorgorjanc Yes, I see the importance. I will propagate the AlphaPeel accuracy tests to AlphaImpute2.