Performs genome-wide marginal association analysis from a GDS file with chunked processing for memory efficiency. Supports standard asymptotic and saddlepoint approximation (SPA) methods.
marginal_scan(
gds_file,
null_model,
sample_id = NULL,
use_SPA = NULL,
chunk_size = 2000L,
mac_cutoff = 1L,
missing_imputation = "mean",
output_csv = NULL,
n_cores = 1L,
verbose = 1
)Character path to GDS file.
A glow_null_model object from fit_null_model().
Optional character vector of sample IDs to override
null_model$sample_id for GDS sample alignment. Default: NULL
(uses null_model$sample_id).
SPA control: NULL (auto: TRUE for binary, FALSE for continuous), TRUE (force SPA, binary only), FALSE (standard only), or "both" (equivalent to TRUE). Default: NULL.
Integer number of variants per processing chunk. Default: 2000L.
Integer minimum minor allele count. Variants with MAC < mac_cutoff are excluded from results. Default: 1L.
Character method for missing genotypes: "mean" or "zero". Default: "mean".
Optional path to write results as CSV. Default: NULL.
Integer (default 1). Number of forked workers over which the
variant chunks are spread (parallel::mclapply). The chunks are
split into n_cores contiguous groups, each worker opens its own
GDS handle and processes its group in order, and the groups are bound in
order, so the result is the same as with one core. Forking is not
available on Windows, where a value above 1 warns and runs serially.
Memory grows with the workers (each holds its own chunk), and a machine
with a threaded BLAS should set OPENBLAS_NUM_THREADS=1 to avoid
oversubscription.
Integer verbosity level (0=silent, 1=progress, 2=detailed).
A data.frame with class "glow_marginal_scan" containing columns: chr, pos, ref, alt, variant_id, MAF, MAC, score, var_score, Z, pvalue. When SPA is active, additional columns Z_SPA and pvalue_SPA are included.
For each chunk of variants, genotype dosages are read from the GDS file using SeqArray, missing values imputed (mean or zero), and marginal score statistics computed from the pre-fitted null model.
Score statistic: \(S_j = G_j^T (Y - \mu_0)\) where \(\mu_0\) are fitted values under the null. Variance is computed using the hat-matrix projection to account for covariates. The Z-statistic is \(Z_j = S_j / \sqrt{Var(S_j)}\).
When SPA is enabled (recommended for binary traits with case-control
imbalance), SPAtest::ScoreTest_SPA() is called per chunk and
SPA-corrected p-values are stored in the pvalue_SPA column.
Complexity: O(n * m) per chunk where n is sample count and m is chunk_size. Memory per chunk is O(n * chunk_size). Total runtime is linear in the number of variants.
Dey, R., Schmidt, E. M., Abecasis, G. R., and Lee, S. (2017). A fast and accurate algorithm to test for binary phenotypes and its application to PheWAS. American Journal of Human Genetics, 101(1), 37-49. doi:10.1016/j.ajhg.2017.05.014
if (FALSE) { # \dontrun{
pheno <- extract_pheno_covar_gds("chr22.gds", covar_names = "sex")
nm <- fit_null_model(pheno$X, pheno$Y, trait = "binary",
sample_id = pheno$sample_id)
results <- marginal_scan("chr22.gds", nm, use_SPA = FALSE, verbose = 1)
head(results)
} # }