[WIP] Parallel linclust

[bbuschkaemper]

WIP

TL;DR

This is a work-in-progress parallelization of the linclust algorithm using a multi-node setup with a clustered filesystem (e.g. GPFS, Lustre) without requiring cross-node communication (e.g. MPI) to address limitations when trying to cluster >1e11 sequences on a single node.

Shared-filesystem parallel linclust

Base problem

To cluster 1e11–1e12 sequences on nodes with <2 TB RAM and <1 PB of scratch we run into an issue; several structures are sized by key space or by the whole database:

structure	at 1e12
seqkey_to_len[dbKeySize], countTable, repSequence (kmermatcher.cpp)	2–4 TB each
resident DBReader::Index[]	~24 TB
assignedCluster[dbSize] (Align2clust.cpp:445)	8 TB
mergeclusters' std::list<size_t>[N] (mergeclusters.cpp:28)	24 TB of empty headers

Plus a time cost: when the k-mer array does not fit, kmermatcher splits and re-extracts every
k-mer per split - roughly 68 full extractions at 1e12.

Basic idea

Replace splitting with one scan plus a shuffle, and partition k-mer space rather than sequence
space. The partition of a k-mer is the low bits of the hashUInt64 score kmermatcher already
computes - a pure function of the k-mer - so every occurrence of a k-mer lands in the same
partition, and a partition can be grouped in isolation. The division is lossless, not
approximate, which is what sequence-space sharding cannot offer.

Coordination is files only. Every worker of a stage runs a byte-identical command line;
identity comes from a fetch_add on a counter file. There is no rank argument and no
node-to-node communication, so a stage maps onto a Slurm array job and workers may join late, die,
or be restarted. Mutual exclusion uses POSIX fcntl(F_SETLKW) whole-file locks, chosen
deliberately over flock, which is node-local on GPFS and Lustre.

Pipeline

Nine new commands driven by one workflow script:

createdbparallel      FASTA -> sequence DB with dense, length-ranked keys
  kmermatcherparallel   extract k-mers into P partition buckets (one wave at a time)
  kmerreduceparallel    per partition: sort, assignGroup -> candidate edges, bucketed by rep key
  alignparallel         per rep-key bucket: merge duplicate copies, align once
  greedycluster         single node, one key-ordered sweep
createrepdb / translatecluster / mergeclusterparallel / translatekeys

Design decisions

1. Length-ranked dense keys. createdbparallel emits sequences longest-first, so
key == global length rank. This removes SORT_BY_LENGTH and its side arrays, lets a companion
index address an entry by key with no resident .index, and makes stock's
longest-first greedy identical to ascending key order. The clustering collapses to a single
left-to-right sweep needing 2 bits per key instead of 8 bytes: 25 GB at 1e11 against 800 GB.

2. The alignment is keyed by representative range, not by k-mer partition. We implemented the
obvious fused design first and measured it: a k-mer partition's pairs are scattered over the whole
key space, so one partition needed 1.61 GB of sequences but read 83.9 GB - 52× amplification,
essentially the whole database per partition. Bucketing surviving edges by representative key
gives each align worker a contiguous slice (measured 1.15×). Two properties fall out for free:
cross-partition duplicate pairs are removed before aligning, and stock's global
per-(pair, diagonal) score accumulation is reproduced exactly rather than approximated.

3. Sizing is derived, not configured. --scratch-budget is the only new scale knob; the
extraction wave count and partition count follow from it and --split-memory-limit. A wave
re-extracts every k-mer but keeps only its slice of partition space, so peak scratch is the whole
shuffle divided by the wave count.

Verification

Against stock run on the same database:

• pass 1 is exact - 0 of 1,000,000 sequences differing;
• end-to-end, 8 of 1,000,000 (0.0008%) sit in a differing cluster, traced to a bounded
  semantic difference in (pair, diagonal) accumulation;
• at 1e9 MGnify sequences, 6,504 of 1e9 (0.00065%), with peak scratch 767.7 GB against a
  1 TB budget;
• output is byte-identical across worker counts, wave counts, thread counts, and after a
  kill-and-resume;
• verified multi-node on JURECA-DC HPC (GPFS), where the multi-node path contributed no divergence of
  its own.

Limitations

• --cov-mode 1 or 2 only. Symmetric coverage modes make linclust select SET_COVER plus
  the count-table rounds; we implement neither, so the command refuses them.
• Protein only. alignparallel rejects nucleotides.
• Parity is measured against createdb --shuffle 0, where it is exact. createdb's default
  32-way shuffle changes which of several equal-length sequences wins a tie, and
  createdbparallel does not reproduce it.
• Output is representative<TAB>member in accessions, not a cluster DB - a per-key index is
  a state no single node can hold at this scale.
• Stock is touched, in 6 files (+150/−77): two defaulted-NULL parameters and guarded
  branches in kmermatcher.{cpp,h}, and parsePrecisionLib de-duplicated into Matcher.cpp.
  Both checked behaviour-preserving - linclust on 1M sequences produces a clusters.tsv
  byte-identical to a reference binary built from the base commit. Separately, Parameters.cpp
  relaxes checkIfDatabaseIsValid so a mkdir race is tolerated when the directory already
  exists; this affects every command with a directory output validator and fixes a real
  Slurm-array race.