Sparse-LD preparation contract

Scope and ownership

gbase::gprep(Glist, task = "sparseld", ...) validates an existing prepared gs_genotypes descriptor, resolves its physical BED rows and marker records, and delegates computation and persistence to gmat. gbase owns the public task, validation, selection, and returned R descriptor. gmat owns the numerical kernel, CSR storage, validation, and LD-score scan. gbits remains the existing packed BED access path. BED remains the genotype store, and neither a dense genotype matrix nor a dense genome-wide LD matrix is constructed.

The input is never prepared implicitly. Unknown tasks and the separate "ldscores" task fail. qgg and sblr are read-only qualification references and are not runtime dependencies.

Sample and marker resolution

Reference samples are selected by explicit ids, then Glist$study_ids, then all Glist$ids. Values must be nonempty, unique, nonmissing, and completely resolved. Caller order is preserved in idsLD, provenance, and the exact zero-based source rows passed to gmat.

Markers are selected by explicit rsids, otherwise all prepared markers. Requested IDs must be unique and completely resolved. The physical matrix is always ordered by file and original BIM order; requested order is recorded in provenance. source_marker_indices_0based supplies the original BED record for every marker, so noncontiguous preparation and LD subsets are safe. LD never crosses a chromosome boundary.

Allele frequencies count BIM field 5 (a1). Frequencies from preparation are reused only when the ordered LD samples equal the preparation sample set. Otherwise the existing gbits packed-count primitive recomputes frequencies for the LD reference rows; an all-missing marker fails. This frequency vector is aligned exactly to the CSR marker order. gmat independently obtains the same reference-sample means while standardizing genotypes.

Numerical contract

For observed allele-1 dosage g, gmat uses \(p=\operatorname{mean}(g)/2\) and \(z=\frac{g-2p}{\sqrt{2p(1-p)}}\). A missing call is standardized zero, equivalent to mean imputation after centering. Monomorphic and all-missing columns have zero scale and become zero; gsuite rejects an all-missing selected reference marker before computation. Pairwise dot products accumulate in double precision and persisted correlations are signed float32 r. A zero denominator produces \(r=0\).

Only strict-upper off-diagonal entries satisfying \(r^2\geq\texttt{r2}\) are stored, in deterministic row/column order; the unit diagonal is implicit. A base-pair window and marker-index window are inclusive and combined with AND when both are positive. Zero disables one window. Both zero are rejected unless allow_full_ld = TRUE because that evaluates every same-chromosome pair.

Ordinary LD score is computed by scanning the validated persisted CSR: \(1+\sum r^2\) over retained off-diagonal neighbors. The implicit diagonal contributes one, thresholded-out entries do not contribute, and no finite- sample correction is applied. Scores follow rsidsLD exactly.

Disk and R descriptors

The canonical gmat format is sparse_ld_csr, format version 1. For prefix P, it owns exactly P.row_ptr.u64.bin, P.col_idx.u32.0based.bin, P.values.f32.bin, P.markers.tsv, and P.meta.txt. The first three contain a strict-upper CSR with zero-based columns and native fixed-width scalar types; the marker sidecar freezes identity and physical source mappings. Exact sizes, row pointers, sorted columns, finite values, dimensions, and marker metadata are validated before publication and on reopen.

Writes use a same-directory staging set. overwrite = FALSE refuses any existing owned target. Explicit overwrite replaces only that five-file set; metadata is published last as the commit marker. Staging and partially published new files are cleaned on failure.

The returned object remains c("gs_genotypes", "Glist"). Its canonical sparseLD member records prefix, format/version, paths, sample/marker counts, nnz, component files, physical one-based rows/cls, aligned AF, marker and sample IDs, zero-based mappings, thresholds/windows, block/thread settings, and provenance. Top-level ldfiles, idsLD, and rsidsLD are the public resource, reference-sample, and ordered-marker identities for that same LD computation; ldscores carries its aligned ordinary scores. Validation requires these views to agree with sparseLD, so they cannot describe a different resource. No external pointer is stored, so saveRDS()/readRDS() is exact and the files can be reopened after reload.

Creation fully validates the published native resource and records checksums and ordinary LD scores. Normal validate_Glist() calls check the serializable descriptor, identity fields, canonical file availability, dimensions, and alignment without hashing or rescanning the complete CSR payload. Use verify_ld() after receiving, copying, or downloading a resource, or when corruption is suspected; it reopens the native resource and verifies its structure and metadata, recorded checksums, and ordinary LD scores. Sparse-LD files are immutable between verification and analysis.

Native consumers open that already-checked identity once per high-level call. The shared gmat reader loads descriptor metadata and row pointers immediately. Regional reads seek directly to selected CSR rows; the first active-column request builds one O(marker count) reverse index that later column windows reuse. Multiple regions, marker sets, traits or chains share the same reader; only the active dense regional matrix or bounded column page is materialized. This keeps repeated full verification out of numerical loops. Consumers do not retain a process-global reader or store an external pointer in Glist, so R serialization and explicit resource ownership remain unchanged.

Independent blocks can be attached after sparse-LD preparation:

blocks <- data.frame(
  block_id = c("1a", "1b"),
  first_marker = c("rs1", "rs101"),
  last_marker = c("rs100", "rs200")
)
Glist <- gbase::define_ld_blocks(Glist, blocks)

The blocks must be ordered, contiguous, exhaustive and confined to one chromosome each. Definition performs one complete native traversal and rejects any retained cross-block correlation. The resulting ldBlocks manifest is tied to the recorded sparse-LD checksums and survives serialization. C++ consumers can then open one session reader and materialize one block at a time; eigensystem truncation and statistical likelihood choices are not stored in Glist.

The core payload layout is readable by the committed sblr CSR reader. The tested qgg 1.1.6 sparse-LD task instead writes fixed-width band files, so compatibility with qgg is numerical rather than a shared persistent format. The familiar msize argument aliases the marker-distance window. With msize and omitted distance/cutoff settings, no base-pair window or r² cutoff is added. Explicit modern settings remain available. ldfiles supplies one output prefix; its optional .ld suffix is removed. getLDscores() returns stored marker-named scores without opening LD files. createLDsets() finds candidate score-based regions; it does not verify independent blocks. No converter or sblr adapter is introduced.

Complexity

Persistent input is the existing \(O(NM/4)\) BED and output is \(O(M + \mathrm{nnz})\). Working memory is bounded by \(O(M+\texttt{source\_N}+\texttt{selected\_N}\,\texttt{block\_size}+\texttt{active\_row\_nnz})\) plus per-row thread-local retained entries. It does not hold \(O(M^2)\) LD or \(O(NM)\) decoded genotypes. Realistic reference LD should be prepared once and reused.