Some clarification questions for benchmarking

[shy218]

I have a few concerns about this research and our goals.

1. What are the packages we want to compare? I currently select the SGKIT, Hail, Glow, Plink, pbdR and SystemML. pbdR and SystemML are two low-level packages, which are not specific to genomics datasets. If we add these two packages, I need to reimplement all the algorithms again from scratch, and the workload is huge. This is not doable in just two months.
2. And also, there is another related question. What is our input format? VCF file or genotype call? If our input format is VCF file, do we need to also write the functions to convert VCF to genotype call? Or could we call the functions from SGKIT? If our input format is genotype call, some packages, such as Plink and Glow cannot be tested. I currently prefer to test on some packages with the same input formats and we select some functions that all the packages have.
3. If I need to implement the algorithms on my own, my implementation may be biased. I think our goal is to test and benchmark the performance of the packages for genomics datasets, not to reimplement and try to optimize the algorithm. I am not an expert in Bioinfo, and maybe some implementation by myself is not optimized.
   Sorry for bothering you time. I am confused after the deep analysis of the packages and algorithms this week.

[hammer]

What are the packages we want to compare?

I think sgkit, Hail, and PLINK are fine to start.

What is our input format?

That's a good question. I think we should use VCF as the starting point and allow for a conversion step for each library into the library's preferred format, with subsequent operations making use of that preferred format.

If I need to implement the algorithms on my own

You should not need to implement any algorithms on your own.

[hammer]

Plink is like a black box to me

If you mean the code is hard to understand, that's fair! In case you mean you don't know where to find the code, it's at https://github.com/chrchang/plink-ng.

[eric-czech]

Hail has both KING and PC-Relate. Should we try to compare to them somehow?

Good point, I think an sgkit.pc_relate comparison to Hail's pc_relate would be excellent. @ravwojdyla has good examples for that in https://github.com/related-sciences/gwas-analysis/blob/rav-pc_relate-tw/notebooks/platform/xarray/pc_relate_perf.py#L3-L7 and perhaps elsewhere?

@eric-czech any thoughts on if/how we should benchmark relatedness estimation and pruning

I don't see much value in doing relatedness pruning since making a representative choice from such small groups is easy by comparison to the same problem with LD. In other words, the family clusters are much easier to isolate than LD clusters and pruning AFAIK typically means just picking someone from a cluster with some ad-hoc routine.

@aktech's great work on pairwise distance computation?

I think the most natural target for this is https://hail.is/docs/0.2/methods/stats.html#hail.methods.row_correlation, which is the same as computing a dense LD matrix for specific inputs (i.e. pearson correlation between variants). That would be a nice one to include since it's probably more applicable to popgen work as well.

[ravwojdyla]

Regarding PC-Relate comparison with Hail, this is a good start: https://discourse.pystatgen.org/t/pc-relate-experiment/51, there are further links in that post.

Also regarding comparison with Hail, as pointed out in related-sciences/gwas-analysis#4 (comment):

it would be interesting to see a perf difference between numpy (and/or Dask Array) and Hail's BlockMatrix.

Hail's BlockMatrix could be compared with Dask's Array, on linear algebra operators. This could be a "low" level benchmark, that builds up to more use case/genetic specific benchmarks.

[shy218]

We can replace "call rate" benchmark to pc-relate, since pc-relate is more interesting to be compared.

[shy218]

I think we should have three axis of benchmarks -- Complexity or Difficulty of Algorithm, Data Factors (Data size) and Computational Factors (scalability). These are consistent with the SLAB paper from Professor Arun Kumar.

[hammer]

I wonder if these benchmarks are overly focused on statistical genetics. @tomwhite are there any unique computational patterns from the population genetics methods you've implemented that might make sense to add to the benchmark?

[tomwhite]

The Garud H stats are the most computationally expensive of the popgen methods that I ported. This is because entries in the samples dimension are hashed and sorted. I don't think there are Hail or PLINK equivalents though, so there wouldn't be anything to compare to.