Design principles for a new GWAS Toolkit

[hammer]

I figured I'd start a topic to collect the design principles that have been implicit in our discussions of what we'd like to build.

Some initial thoughts:

• Written in Python
• Built with and conforming to the code style and extension APIs of popular libraries from the PyData ecosystem such as NumPy, Pandas, Dask, Xarray, Zarr, and others.
• In particular, we should not have to build our own data frame or n-dimensional array APIs. We should inherit those APIs from upstream libraries.
• Make use of type checking as much as possible.
  • Cf Support static type analysis  pydata/xarray#3967, Type hinting / annotation (PEP 484) for ndarray, dtype, and ufunc numpy/numpy#7370
• Able to make use of CPU instruction set extensions like AVX
• Able to use GPUs.
• Able to scale up to multiple nodes.
• Have a uniform API that can dispatch to multiple backends depending on the workload.
• Make use of standardized container file formats like Zarr.
• Work well on the public cloud.

[hammer]

Some questions to answer

• What organisms should we target?
• Should we try to accommodate single cell workloads too?

[hammer]

Some optimizations that would be nice

• Able to operate on sub-byte data representations without requiring a custom API.
• Able to optimize over a pipeline of operations (noted as important in A comparative evaluation of systems for scalable linear algebra-based analytics (2018)).

[hammer]

(Post by @jeromekelleher)

What organisms should we target?

I would say "all of them". Anything that we specialise for haploids, diploids, haplodiploids, etc will end up being limiting and short-sighted. This should be a general-purpose genomics library (that works well in the common diploid case).

Should we try to accommodate single cell workloads too?

I would say "no" as it seems like mission creep, but then I don't do any single cell so I don't know.

[hammer]

(Post by @alimanfoo)

Able to operate on sub-byte data representations without requiring a custom API.

FWIW I previously tried to avoid anything sub-byte in scikit-allel, mainly because numpy is byte-atomic. Also bit-shifting is not something I find easy to code, I know it's bread and butter to lots of folks coming from a more computational background, but can be less natural to hackers coming from a more biological background (like me :-)). Also I found that whole-byte representations with compression were as compact or more so than sub-byte representations, so using sub-byte representations just to save space doesn't always make sense. But these are small points, I don't want to shut any ideas down here, just some thoughts.

[hammer]

A discussion on Dask configuration happening at Related Sciences right now between @eczech and @ravwojdyla has reminded me of another design principle

• Good default configuration parameters, with limited manual tuning required for new backends or workloads.

As noted in https://discourse.pystatgen.org/t/a-comparative-evaluation-of-systems-for-scalable-linear-algebra-based-analytics-2018/15, bad default configuration parameters and lots of workload-specific tuning are big drawbacks of Spark.

In particular, we should have default configuration parameters that work well on a single node, as @ravwojdyla has started to collect at related-sciences/gwas-analysis#27.

[hammer]

(Post by @ravwojdyla)

If I may add some thoughts here, from working on some other processing frameworks, things that I find useful:

• easy/common problems should be easy to solve (build in support), hard/esoteric problems should be possible to solve (using build in blocks)
• as much determinism as practically feasible, users should not be surprised about the frameworks behavior and execution plan/steps (I have seen and been victim to this for example in Hail)
• it should be easy to test your pipeline, so having a scaffolding for tests, and maybe data generators should enable people to write tests and maintain reasonable trust in their codebase (this has been a problem in Spark)
• reproducibility - for example a failure in distributed mode should be possible to reproduce in "local" mode (Apache Beam has this problem)
• observability - there is nothing worse then a pipeline that does not make progress and you can't easily check whether you should give up, whether there is an issue (struggler) etc. This usually has been solved at scale with counters, gauges, dashboards
• debugability - in case of bugs and issues, in perfect scenario, in local mode, I should be able to set breakpoints and step through the execution
• there should be an efficient local mode (Apache Beam has this problem too)
• consistent and typed API, I know this thing can be hard to achieve in Python. The benefits of typing are unquestionable. I think it's good to start typing early, before there is too much code that is "impossible" to type in python. Would be great if at least the public API was typed.
• thinking about the execution time of some of these pipelines, I wonder if it's also important to think about some support for checkpointing and restarts
• non cryptic serde failures would be nice to have