Requirements for UKB GWAS

[eric-czech]

To run a basic GWAS on UKB data, here are some of the operations we'll need support for:

• bgen reader (BGEN reader implementation using bgen_reader sgkit-bgen#1)
• plink reader (Pysnptools reader implementation sgkit-plink#1, Suboptimal parallelism sgkit-plink#6)
• Variant allele frequency/count (https://github.com/pystatgen/sgkit/issues/29)
• Variant call rate/count (https://github.com/pystatgen/sgkit/issues/29)
• Variant HWE test (https://github.com/pystatgen/sgkit/issues/28)
• Sample call rate/count (https://github.com/pystatgen/sgkit/issues/29)
• An is_autosome function to filter variants by
• A function to convert genotype probabilities to hard calls (https://github.com/pystatgen/sgkit/issues/346)
• A linear regression function (https://github.com/pystatgen/sgkit/pull/52)
• A variant annotation function like vep. There are plenty of other ways to get this but an internal function would be great.
• A phenotype normalization pipeline. I don't expect much of this to become part of sgkit, but there might be some generalizable phenotype-specific functions that are worth considering for inclusion.

There may be a few more beyond that, but I think anything remaining should be reasonable with Xarray/Dask alone.

[tomwhite]

I'd be happy to work on variant/sample stats (#29) if no one else is working on them.

[hammer]

@eric-czech how are you thinking about LD estimation/pruning, population structure estimation/pruning, and relatedness estimation/pruning? Does REGENIE include in the implementation some means of estimating these things as covariates for the regression, or are you just thinking of those operations as optimizations that can be implemented later?

[hammer]

To answer my own question, there are 3 stages to our work with UK Biobank

• Stage 1: per-variant linear regression with provided population structure and kinship estimates. No LD pruning needed.
• Stage 2: whole-genome regression with provided population structure and kinship estimates. LD pruning needed.
• Stage 3: do our own population structure and kinship estimation.

[hammer]

A phenotype normalization pipeline.

I've been thinking about this one too. I think we're going to feel Dask's poor handling of nested data when working with phenotypes, and I'd prefer to keep Spark out of this project as a dependency, so I think we put that code into a separate repo if we find we do need Spark.

[hammer]

A variant annotation function like vep. There are plenty of other ways to get this but an internal function would be great.

File an issue to track?

[eric-czech]

File an issue to track?

https://github.com/pystatgen/sgkit/issues/112

[hammer]

A phenotype normalization pipeline.

@eric-czech I just had a nice chat with @zietzm and @ntatonetti who are at Columbia and are experts in handling complex phenotypes and running many GWAS against them.

They're interested in using sgkit and possibly contributing back, particularly on the phenotype side.

Would you be open to making https://github.com/related-sciences/ukb-gwas-pipeline-nealelab public soon and potentially working with @zietzm to factor the phenotype handling code into its own repo, maybe something like sgkit-pheno or phenokit?

[eric-czech]

Would you be open to making https://github.com/related-sciences/ukb-gwas-pipeline-nealelab public soon and potentially working with @zietzm to factor the phenotype handling code into its own repo, maybe something like sgkit-pheno or phenokit?

For sure! Looking forward to seeing how we can better integrate phenotypes.

[eric-czech]

FYI @zietzm / @ntatonetti (cc: @hammer) the phenotype prep code we're currently using (via PHESANT) is here: ukb-gwas-pipeline-nealelab#phenotype_prep.smk.

There is little to it yet other than running some messy, very inefficient R code to produce ~75 phenotypes that I wanted to attempt to validate against first. It would be great to hear your thoughts on how we might better define these as well as improve the mechanics of how we're creating them. I'm particularly interested in ICD code management since this pipeline doesn't address that.