Without the debug swith gemma runs at the same speed with 128% CPU. That won't help much.

Compiling with optimization can be low hanging fruit - despite the fact that we seem to be IO bound at 128% CPU. Still, aggressive compiler optimizations may make a difference. The current build reads:

GEMMA Version    = 0.98.6 (2022-08-05)
Build profile    = /gnu/store/8rvid272yb53bgascf5c468z0jhsyflj-profile
GCC version      = 14.3.0
GSL Version      = 2.8
OpenBlas         = OpenBLAS 0.3.30  - OpenBLAS 0.3.30 DYNAMIC_ARCH NO_AFFINITY Cooperlake MAX_THREADS=128
arch           = Cooperlake
threads        = 96
parallel type  = threaded

this uses the gemma-gn2 package in

• https://git.genenetwork.org/guix-bioinformatics/tree/gn/packages/gemma.scm#n27

which is currently not built with arch optimizations (even though Cooperlake suggests differently). Another potential optimization is to use a fast malloc library. We do, however, already compile with a recent gcc, thanks to Guix. No need to improve on that.

Rather than focussing on gzip, another potential improvement is to use lmdb with mmap. We am not going to upgrade the original gemma code (which is in maintenance mode). We are going to upgrade the new pangemma project instead:

• https://git.genenetwork.org/pangemma/

Reason being that this is our experimental project.

So I just managed to build pangemma/gemma in Guix. Next step is to introduce lmdb genotypes. Genotypes come essentially as a matrix of markers x individuals. In the case of GN geno files and BIMBAM files they are simply stored as tab delimited values and/or probabilities. This happens in

src/param.cpp
1261:void PARAM::ReadGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) {
1280:void PARAM::ReadGenotypes(vector<vector<unsigned char>> &Xt, gsl_matrix *K,

calling into

gemma_io.cpp
644:bool ReadFile_geno(const string &file_geno, const set<string> &setSnps,
1752:bool ReadFile_geno(const string file_geno, vector<int> &indicator_idv,
1857:bool ReadFile_geno(const string &file_geno, vector<int> &indicator_idv,

which are called from gemma.cpp. Also lmm.cpp reads the geno file in the AnalyzeBimbam function (see file_geno):

src/lmm.cpp
61:  file_geno = cPar.file_geno;
1664:  debug_msg(file_geno);
1665:  auto infilen = file_geno.c_str();
2291:    cout << "error reading genotype file:" << file_geno << endl;

Note that also SNPs are read from a file (see file_snps). We already have an lmdb version for that!

So, reading genotypes happens in multiple places. In fact, it is read 1x for computing K and 2x for GWA. And it is worth than this because LOCO runs GWA 20x rereading the same files. Reading it once using lmdb should speed things up.

We'll start with the 30G 143samples.percentile.bimbam.bimbam-reduced2 file. To convert this file into lmdb we only do this once. We want to track both column and row names in the same lmdb and we will use a meta JSON record for that. On the command line we'll state wether the genotypes are stored as char or int. Floats will be packed into either of those. We'll expirement a bit to see what the default should be. A genotype is usually a number/character or a probability. In the latter case we don't have to have high precison and can choose to store an index into a range of values. We can also opt for Float16 or something more ad hoc because we don't have to store the exponent.

But let's start with a standard float here, to keep things simple. To write the first version of code I'll use a byte conversion:

./bin/geno2mdb.rb BXD.geno.bimbam --eval '{"0"=>0,"1"=>1,"2"=>2,"NA"=>-1}' --pack 'C*' --geno-json BXD.geno.json

The lmdb file contains a metadata record that looks like:

{
  "type": "gemma-geno",
  "version": 1,
  "eval": "G0-2",
  "key-format": "string",
  "rec-format": "C*",
  "geno": {
    "type": "gn-geno-to-gemma",
    "genofile": "BXD.geno",
    "samples": [
      "BXD1",
      "BXD2",
      "BXD5",
etc.

i.e. it is a self-contained, efficient, genotype format. There is also another trick, we can use Plink-style compression with

./bin/geno2mdb.rb BXD.geno.bimbam --eval '{"0"=>0,"1"=>1,"2"=>2,"NA"=>4}' --geno-json BXD.geno.json --gpack 'l.each_slice(4).map { |slice| slice.map.with_index.sum {|val,i| val << (i*2) } }.pack("C*")'

reducing the original uncompressed BIMBAM from 9.9Mb to 2.7Mb. This is still a lot larger than the gzip compressed BIMBAM, but as I pointed out earlier the uncompressed version is faster by a wide margin. Compressing the lmdb file gets it in range of the compressed BIMBAM btw. So that is always an option.

Next we create a floating point version. That reduces the file to 30% with

geno2mdb.rb fp.bimbam --geval 'g.to_f' --pack 'F*' --geno-json bxd_inds.list.json

and if we compress the probabilities into a byte reduces the file to 10%:

geno2mdb.rb fp.bimbam --geval '(g.to_f*255.0).to_i' --pack 'C*' --geno-json bxd_inds.list.json

And now the compressed version is also 4x smaller. We'll have to run gemma at scale to see what the impact is, but an uncompressed 10x reduction schould have an impact on the IO bottle neck. Note how easy it is to try these things with my little Ruby script.

• https://github.com/genetics-statistics/gemma-wrapper/blob/master/bin/geno2mdb.rb

Rather than writing new code in C++ I proceeded embedding guile in pangemma. If it turns out to be a performance problem we can always fall back to C. Here we show a simple test witten in guile that gets called from main.cpp:

• https://git.genenetwork.org/pangemma/commit/?id=5b6b5e2ad97b4733125c0845cfae007e8094a687

GEMMA::BatchRun reads files and executes (b gemma.cpp:1657) cPar.ReadFiles() ReadFile_anno ReadFile_pheno ReadFile_geno (gemma_io.cpp:652) - first read to fetch SNPs info, num (ns_tset) and total SNPs (ns_total). - it also does some checks Note: These can all be handled by the lmdb files. So it saves one run.

Summary of Mutated Outputs:

• indicator_snp: Binary indicators for which SNPs passed filtering
• snpInfo: Complete metadata for all SNPs in the file
• ns_test: Count of SNPs passing filters

checkpoint("read-geno-file",file_geno);

Next start LMM9 gemma.cpp:2571 ReadFile_kin EigenDecomp_Zeroed 2713 CalcUtX(U, W, UtW); 2714 CalcUtX(U, Y, UtY); CalcLambda CalcLmmVgVeBeta CalcPve cPar.PrintSummary() debug_msg("fit LMM (one phenotype)"); cLmm.AnalyzeBimbam lmm.cpp:1665 and LMM::Analyze lmm.cpp:1704

Based on LLM code analysis, here's what gets mutated in the 'LMM' and Param class:

From above it should be clear that, if we have the genotypes and snp annotations in lmdb, we can skip reading the genotype file the first time. We can also rewrite the 'analyze' functions to fetch this information on the fly.

Note that OpenBLAS will have to run single threaded when introducing SNP-based threads.

From above it can be concluded that we can batch process SNPs in parallel. The only output is sumStat and that is written at once at the end. So, if we can collect the sumStat data without collision it should just work.

Interestingly both Guile and C++ have recently introduced fibers. Boost.Fiber looks pretty clean:

#include <boost/fiber/all.hpp>
#include <vector>
#include <iostream>

namespace fibers = boost::fibers;

// Worker fiber
void compute_worker(int start, int end,
                   fibers::buffered_channel<int>& channel) {
    for (int i = start; i < end; ++i) {
        channel.push(i * i);
    }
}

void parallel_compute_fibers() {
    fibers::buffered_channel<int> channel(100);

    // Spawn fibers
    fibers::fiber f1([&]() {
        compute_worker(0, 100, channel);
        channel.close();  // Signal completion
    });

    fibers::fiber f2([&]() {
        compute_worker(100, 200, channel);
    });

    // Collect results
    std::vector<int> results;
    int value;
    while (fibers::channel_op_status::success == channel.pop(value)) {
        results.push_back(value);
    }

    f1.join();
    f2.join();

    std::cout << "Total results: " << results.size() << std::endl;
}

Compare that with guile:

(use-modules (fibers)
             (fibers channels))

;; Worker that streams individual results
(define (compute-worker-streaming start end result-channel)
  (let loop ((i start))
    (when (< i end)
      (put-message result-channel (* i i))
      (loop (+ i 1))))
  ;; Send completion signal
  (put-message result-channel 'done))

;; Collector fiber
(define (result-collector result-channel num-workers)
  (let loop ((results '())
             (done-count 0))
    (if (= done-count num-workers)
        (reverse results)
        (let ((msg (get-message result-channel)))
          (if (eq? msg 'done)
              (loop results (+ done-count 1))
              (loop (cons msg results) done-count))))))

(define (parallel-compute-streaming)
  (run-fibers
   (lambda ()
     (let ((result-channel (make-channel)))

       ;; Spawn workers
       (spawn-fiber
        (lambda () (compute-worker-streaming 0 100 result-channel)))
       (spawn-fiber
        (lambda () (compute-worker-streaming 100 200 result-channel)))

       ;; Collect results
       (result-collector result-channel 2)))))

The Boost fiber is a relatively mature library now, with about 8+ years of development and real-world usage. Interestingly Boost.fibers has work stealing built in. We can look at that later:

• https://www.boost.org/doc/libs/1_66_0/libs/fiber/doc/html/fiber/worker.html

What about LOCO? Actually we can use the same fiber strategy for each chromosome as a per CHR process. We can set the number of threads differently based on chromosome SNP num, so all chromosomes take (about) the same time. Later, we can bring LOCO into one process with the advantage that the genotype data is only read once. In both cases the kinship matrices are in RAM anyway.