Genotype call array to dosage

[alimanfoo]

Raising this issue to discuss implementing a function to convert a genotype call array to a dosage array.

In general, the output array should have at least two dimensions, with the first two dimensions being (variants, samples). The array elements give the dosage of each allele, i.e., how many copies of an allele are carried by the individual.

Some questions for discussion

• How do we handle biallelic and multiallelic variants?
• How do we handle missing genotype calls and avoid creating a bias towards a particular allele (e.g., the reference allele)?
• What should the output dtype be? (int or float or either)? (Some programs think of dosage as a continuous variable.)

Some related functions (not necessarily what we want to copy, but for reference):

• scikit-allel has to_n_ref, to_n_alt, to_allele_counts
• skallel prototype has genotypes_3d_to_allele_counts, genotypes_3d_to_allele_counts_melt, genotypes_3d_to_major_allele_counts

[alimanfoo]

Perhaps this issue should be merged with #3?

[tomwhite]

In #38 I have created an implementation that

• Is biallelic only
• Represents missing calls with NaN in call/dosage and a corresponding mask in call/dosage_mask
• Uses float types

[hammer]

Is biallelic only

I don't want to complicate things too early but as WES/WGS becomes more common than genotyping arrays, multiallelic sites become more common.

Should we strive for simplicity early (i.e. biallelic only) so that we can solve for other issues in library design now, or should we design for complexity up front (i.e. multiallelic sites) so that we don't have to revise things later?

I can see arguments for both approaches, so I'm happy to defer to the team's desires. I think @alimanfoo already works primarily with WGS data, so getting his insights would be valuable.

[jeromekelleher]

In general I think we need to assume everything is multiallelic, from the start. However, is it clear that we need to worry about this for dosages? This is a stupid question, but can we get dosage data for anything other than genotype arrays?

[eric-czech]

can we get dosage data for anything other than genotype arrays?

It's pretty common in GWAS to compute dosage based on genotype likelihood from WES/WGS for regression rather than hard calls, so that's one example at least.

[jeromekelleher]

OK, seems to me like we have to grasp the multi-allelic nettle from the start then.

[alimanfoo]

+1 for thinking multiallelic from the start, will make some entomologists very happy :-)

(And also +1 for thinking about any ploidy from the start - will hopefully make us some botanist friends :-))

[alimanfoo]

can we get dosage data for anything other than genotype arrays?

It's pretty common in GWAS to compute dosage based on genotype likelihood from WES/WGS for regression rather than hard calls, so that's one example at least.

I'm not an expert here, but I was wondering if we should consider separately these two different cases:

• genotype call array -> genotype (per-call) allele counts array
• genotype likelihoods array -> genotype dosage array

Genotype dosage is in general I believe modelled as a continuous variable, i.e., should have a float data type and can take any value between 0 and 2 (for diploid biallelic). I'm guessing genotype dosage is rarely calculated for multiallelic variants.

Perhaps this issue was badly named and/or should be split into two?

[alimanfoo]

Just to add, perhaps there are in fact three cases:

• genotype call array -> genotype (per-call) allele counts array
• genotype likelihoods array (from WES/WGS) -> genotype dosage array
• genotype probabilities array (from genotyping chip + imputation) -> genotype dosage array

The last two cases may be structurally identical in terms of how we represent the data, but are conceptually slightly different and wanted to check these are all paths we want to cover?

[alimanfoo]

Also mentioning one other consideration, handling of missing data. This might not be the best place to discuss, but I'm conscious that converting genotype calls to some kind of allele dosage-like representation can induce a bias in the presence of missing data.

This may not ever come up when dealing with data from chip + imputation because there are no missing calls. But will come up when dealing with WES/WGS data.

The obvious thing to avoid is creating a bias towards the reference allele. E.g., if you assume biallelic sites, and you compute dosage of the alternate allele, then both a hom ref call and a missing call end up with the same dosage (0).