`gs_genotypes` preparation contract

Bounded scope and ownership

gbase::gprep() implements the bounded "prepare", "sparseld", and "eigen" tasks. For preparation, gbase owns argument validation, orchestration, and the returned R descriptor; gmat owns strict BIM/FAM parsing and exact ordered multi-file identity; gbits owns BED header, size, offset, packed-code, and threaded count logic. The existing SNP-major PLINK BED files remain the sole persistent genotype representation. qgg is a read-only behavioural and scientific reference and is not a runtime dependency.

Sparse-LD orchestration is specified separately in sparse-LD preparation. gbase::gprep() does not add separate LD-score tasks: gbase::ld_scores(), ld_region(), and ld_eigen() consume the augmented Glist directly. Genetic-map import, GRMs, association, scoring, and other qgg tasks are not implemented by gbase::gprep(). Marker filtering is provided separately by gbase::gfilter(). Unknown preparation tasks fail explicitly.

Files and identities

Every BED/BIM/FAM vector is nonempty and has identical length. Omitted companions are derived by replacing only a final case-insensitive .bed suffix. Every path is normalized and must name an existing regular file. gbits requires the exact header 6c 1b 01 and exact size $3+M\lceil N/4\rceil$; individual-major, truncated, and trailing-byte inputs fail.

The first FAM supplies ids and the complete physical sample order. Individual IDs are nonempty and globally unique. gmat requires every component to have the same exact ordered (family_id, individual_id) sequence; no set-only check or positional permutation is used. With ids = NULL, all samples contribute and study_ids is NULL. A supplied vector is nonempty, unique, nonmissing, fully resolved, retained in caller order, and used only as the summary sample set; ids remains complete.

Every BIM row has six fields, a nonempty chromosome and marker ID, finite nonnegative map position, positive physical position, and distinct nonempty alleles. Source marker IDs are globally unique across every supplied BIM, including unselected rows. BIM order is canonical. A supplied rsids vector is validated completely, then filtered in file/BIM order; caller order is intentionally not retained because the public descriptor remains a file-partitioned physical view. source_marker_indices_0based records every original BIM/BED index, including noncontiguous selections. BED is never rewritten.

Alleles and summaries

a1 is BIM field 5 and a2 BIM field 6. Counts and frequencies always refer to copies of a1:

PLINK code	Copies of `a1`	Count
`00`	2	`n2`
`10`	1	`n1`
`11`	0	`n0`
`01`	missing	`nmiss`

For $N_{\mathrm{obs}}=n_0+n_1+n_2$, the returned values are $\mathrm{af}=\frac{2n_2+n_1}{2N_{\mathrm{obs}}}$, $\mathrm{af}_1=\mathrm{af}$, $\mathrm{af}_2=1-\mathrm{af}$, $\mathrm{maf}=\min(\mathrm{af},1-\mathrm{af})$, $\mathrm{hom}=(n_0+n_2)/N_{\mathrm{obs}}$, and $\mathrm{het}=n_1/N_{\mathrm{obs}}$. All-missing selected markers fail with their marker ID. Empirical frequencies are not clipped: monomorphic values are exactly zero or one. This deliberately differs from the tested qgg 1.1.6 reference, whose summaryBED() clips frequencies to [0.00001, 0.99999] and substitutes 0.5 for all-missing markers.

Object schema and persistence

The result is a plain list with class c("gs_genotypes", "Glist"). It carries study, normalized bedfiles/bimfiles/famfiles, complete ids, optional study_ids, n, nchr, m, chromosome/file-wise rsids, mchr, chr, pos, map, a1, a2, cpra, af, af1, af2, maf, nmiss, hom, het, n0, n1, n2, assembly = NULL, physical sample/marker mappings, format and ordering labels, and provenance. It contains neither dense genotypes nor a mandatory external pointer and therefore survives an exact saveRDS()/readRDS() round trip.

gsim has a native bounded BED reader for Glist-backed simulation. It resolves selected sample and marker IDs against FAM/BIM files, checks allele orientation, and reads selected BED blocks. A custom deterministic getG_fun remains an explicit alternative. Check the installed gsim version before use.

Chromosome marker sets

Use gbase::getChromosomeSets() to extract chromosomes directly from Glist metadata, without reading BED or LD:

sets <- getChromosomeSets(Glist)

Each chromosome supplies one named, disjoint set covering its markers. Chromosomes retain their first-occurrence order and members retain reference/BIM order. Grouping uses chromosome labels rather than BED file boundaries, so it also works with a single file containing all chromosomes. Pass these sets to glma(..., sets=sets) for chromosome variance components.

Summary-statistics checking and imputation

Keep gaps when extracting database statistics, then prepare all trait columns against the eligible markers in Glist:

stat <- gact::getMarkerStat(GAlist, studyID=c("GWAS1", "GWAS2"), rm.na=FALSE)
stat <- checkStat(Glist, stat, ldcheck=TRUE, impute=TRUE, nthreads=4)
attr(stat, "ldcheck")$original
attr(stat, "imputation")$r2

The same call accepts a summary data frame or a numeric marker-by-trait matrix with marker row names. For a data frame, Z is taken from stat/z, derived from b/seb, or selected with columns=. Alleles are harmonized once; matrix/list alleles can be supplied in the marker_alleles attribute. Without allele metadata, scores must already use Glist’s a1 coding. Ordinary data-frame checkStat() calls retain the established behavior when both options are off. Use impute(Glist, stat) for prediction alone on already aligned input.

Accepted observed Z-scores are preserved. LD-rejected observations are retained in the diagnostic table, masked, and then eligible for imputation. Reference QC exclusions stay excluded; predictions below r2=0.6 remain missing. Diagnostics are attributes on the ordinary input format. Known sample sizes are retained; missing sizes and effect/SE/P fields for new or repaired scores are not inferred. Analyses using imputed scores must account for prediction uncertainty and supply any required sample-size or effect information.

For several ancestries, label samples on the existing descriptor:

Glist <- setReferenceGroups(Glist, list(EUR=eur_ids, EAS=eas_ids))
stat <- checkStat(Glist, stat, ldcheck=TRUE, impute=TRUE)

Groups are disjoint subsets of the prepared samples and are stored as one-based integer indices. Gact carries study ancestry in attr(stat, "ancestry"); an explicit ancestry=c(GWAS1="EUR", GWAS2="EAS") overrides it. Several reference groups require one ancestry choice per trait. A preselected, unlabeled single reference remains supported for one common ancestry.

An explicit ancestry_weights=c(EUR=0.8, EAS=0.2) uses Rmix = 0.8 R_EUR + 0.2 R_EAS; a group-by-trait matrix allows different weights per column. Weights are squared standardized coefficients for independent meta-analysis Z-scores, constant across markers. They are not inferred from the scores and are not a general model of individual admixture or changing study composition/overlap. Allele-frequency differences are checked for pure ancestry columns; covariance weights do not identify a pooled meta-analysis frequency. Reference MAF, missingness and HWE are checked within each contributing ancestry, without pooling populations for HWE.

Imputation uses ridge-conditioned complete regional BED correlation, with lambda=0.1, 1 Mb cores and 250 kb buffers. The default max_markers=2048 bounds local decoded genotypes and dense matrices; dense cores split and their buffers shorten. Matching ancestry and missingness masks reuse factors across traits. No imputed prediction becomes a predictor, no genome-wide LD matrix is built, and no LD or eigenvector files are saved. Stored sparse-LD thresholds are not used. Reference sample selections and allele coding must match the intended downstream LD model.

LD checking uses deterministic repeated held-out partitions, a truncated eigensystem and a default threshold of 5e-8. Near-perfect LD clusters stay in one partition so their members cannot predict each other. This is a DENTIST-style adaptation with bounded windows and clustered partitions; it is not a reproduction of the original program. All checking finishes before imputation, using a second bounded reference pass where predictions are needed. P values and reference prediction quality are approximate. The original DENTIST study recommends at least 5,000 reference individuals. Small deterministic tests verify the computation; genome-scale speed and mismatch calibration have not been established for this new route.

The formulas follow ImpG-Summary and the consistency checks follow DENTIST.

Regional eigenvalue preparation

regionSets <- list(region1 = Glist$rsids[[1]][1:100])
Glist <- gprep(Glist, task = "eigen", sets = regionSets, nthreads = 4)
Glist$eigen$values$region1
saveRDS(Glist, "Glist-with-spectra.rds")

sets is a nonempty uniquely named list of marker-ID vectors. All markers must exist in Glist, duplicates are rejected, and each set must belong to one chromosome. Sets may overlap and need not be contiguous. Marker membership is stored in reference/BIM order. The samples default to study_ids when present, otherwise all ids; an explicit ids argument selects and records another reference sample set. At least two samples are required.

gmat reads unit-norm centred regional dosage columns through gbits. Missing calls are imputed to the observed marker mean; constant or all-missing markers fail with their ID. gsolve computes dense products and eigenvalues only. The spectrum is that of full within-set correlation, without sparse-LD thresholds or distance windows. No eigenvectors, clipping, statistical test, or persistent LD files are produced.

For N selected samples and M markers, the smaller Gram matrix is used. When M exceeds N, the sample-space spectrum is padded with M-N exact zeros and sorted to retain the full M-value spectrum. Floating-point null eigenvalues can differ from a full marker-space decomposition at rounding level. The default eigen_control$maximum_dimension = 2048 applies to min(N, M), not M. One region is processed at a time with O(NM + min(N,M)^2) scratch storage; only O(M) values and region identities remain on Glist. nthreads controls gbits decoding, while the current portable gsolve product and Eigen spectrum backends are serial and do not require OpenBLAS. No global thread setting is changed.

Glist$eigen records named sets, descending values, native diagnostics, selected sample IDs and offsets, controls, and provenance. It survives ordinary RDS persistence. validate_Glist() checks its identities and spectrum shape; it does not recompute the eigenvalues. If a downstream analysis removes markers or changes the reference samples, prepare matching spectra again. glma::prepare_vegas(Glist) reuses the stored spectra; sets may select named regions but must preserve their full membership. Its prepared$markers holds the union of requested markers in reference order. glma::vegas() selects checked statistic rows by marker ID, so the data may contain additional markers or another row order. All required markers and their reference alleles must be present; no QC or allele harmonization is repeated during fitting. Constant positive weights within each set are supported. For nonuniform weights, use glma’s existing sparse-LD preparation on a Glist without the unweighted eigen resource. Save the prepared object to evaluate subsequent traits without reopening the original BED files. gsea re-exports both functions.

Complexity and execution

gbits scans packed marker records directly. Selected original offsets are checked with 64-bit arithmetic. Threaded execution assigns deterministic static marker partitions and never calls R from workers. Temporary native storage is $O(T\lceil N/4\rceil)$ packed bytes for T workers plus $O(M)$ counts and selection mappings; preparation does not build an $O(NM)$ decoded matrix. The separate regional eigen task uses the bounded regional scratch described above.

gs_genotypes preparation contract