Scalable Generalized Mixed Models in PheWAS using SAIGEgds
Sep, 2026
1 Introduction
A phenome-wide association study (PheWAS) is a powerful tool for the discovery and replication of genetic associations across many phenotypes. Large biobanks provide deep genomic and phenotypic data, but testing thousands of traits against whole-genome variants remains computationally demanding. SAIGE (scalable and accurate implementation of generalized mixed models [2]) addresses case-control imbalance and sample relatedness in association studies.
SAIGEgds is a high-performance implementation of the SAIGE framework using Genomic Data Structure (GDS) files [1,3]. Its optimized C++ code takes advantage of sparse genotype dosages and supports integer genotypes and numeric imputed dosages. SAIGEgds supports binary, quantitative and time-to-event outcomes, as well as single-variant and set-based association tests. This vignette focuses on the binary-trait, single-variant workflow; time-to-event outcomes are covered in the companion survival vignette.
Benchmarks using UK Biobank White British genotype data (\(N=430\mathrm{K}\)) with coronary heart disease and simulated case-control outcomes showed that SAIGEgds was 5 to 6 times faster than the SAIGE R package for null-model fitting and p-value calculation. Together with high-performance computing (HPC) clusters or cloud resources, it provides an efficient pipeline for biobank-scale PheWAS.
2 Installation
- Bioconductor repository
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("SAIGEgds")The BiocManager::install() approach may build the package from source. In that case, make and a suitable compiler toolchain must be installed; see the R FAQ for platform-specific instructions.
3 Workflow
library(SeqArray)
library(SAIGEgds)
# the genotype file in SeqArray GDS format (1000 Genomes Phase1, chromosome 22)
(geno_fn <- seqExampleFileName("KG_Phase1"))## [1] "/Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/library/SeqArray/extdata/1KG_phase1_release_v3_chr22.gds"3.1 Preparing SNP data for the genetic relationship matrix
# open a SeqArray file in the package (1000 Genomes Phase1, chromosome 22)
gds <- seqOpen(geno_fn)The LD pruning is provided by snpgdsLDpruning() in the package SNPRelate:
library(SNPRelate)## SNPRelate -- supported by Streaming SIMD Extensions 2 (SSE2)set.seed(1000)
snpset <- snpgdsLDpruning(gds)## SNV pruning based on LD:
## Calculating allele counts/frequencies (19773 variants) ...
##
[..................................................] 0%, ETC: --- (1/1)
[==================================================] 100%, used 1s (1/1)
[==================================================] 100%, complete, 1s
## # of selected variants: 9,043
## Excluding 10,730 SNVs (monomorphic: TRUE, MAF: 0.005, missing rate: 0.01)
## # of samples: 1,092
## # of SNVs: 9,043
## using 1 thread/core
## sliding window: 500,000 basepairs, Inf SNPs
## |LD| threshold: 0.2
## method: composite
## Chrom 22: |====================|====================|
## 14.58%, 2,883 / 19,773 (Thu Sep 24 20:16:22 2026)
## 2,883 markers are selected in total.str(snpset)## List of 1
## $ chr22: int [1:2883] 1 8 9 10 12 16 17 18 20 23 ...snpset.id <- unlist(snpset, use.names=FALSE) # get the variant IDs of a LD-pruned set
head(snpset.id)## [1] 1 8 9 10 12 16Create a genotype file for the genetic relationship matrix (GRM) using the LD-pruned SNP set:
grm_fn <- "grm_geno.gds"
seqSetFilter(gds, variant.id=snpset.id)## # of selected variants: 2,883# export to a GDS genotype file without annotation data
seqExport(gds, grm_fn, info.var=character(), fmt.var=character(), samp.var=character())## ##< 2026-09-24 20:16:22
## Export to 'grm_geno.gds':
## sample.id (1,092) [md5: bd2a93b49750ae227793ed23c575b620]
## variant.id (2,883) [md5: 37e34ffd0328a6bb5bcedb53dddad2e3]
## position [md5: 9f2c5f81ad3709aa4e8c4549a82d2175]
## chromosome [md5: 7b4b80eaf4724d87d0d9b73aab237580]
## allele [md5: 20bd0aa343398d59112bd4b58f562f68]
## genotype [md5: be854ad6089ef332349ff25a99674949]
## phase [md5: 94c1e11249a7e2a4426b02bf852af5cb]
## annotation/id [md5: 2020922f9380eddb5e5f8e271ade54e2]
## annotation/qual [md5: 5d14e0a74e09172e9192c42b646c1559]
## annotation/filter [md5: 9bcb8f5aa38515d50d47129d3c57a542]
## Done. # 2026-09-24 20:16:27
## Optimize the access efficiency ...
## Clean up the fragments of GDS file:
## open the file 'grm_geno.gds' (343.6K)
## # of fragments: 101
## save to 'grm_geno.gds.tmp'
## rename 'grm_geno.gds.tmp' (343.2K, reduced: 444B)
## # of fragments: 64
## ##> 2026-09-24 20:16:27If genotypes are split by chromosomes, seqMerge() in the SeqArray package can be used to combine the LD-pruned SNP sets.
# close the file
seqClose(gds)3.2 Fitting the null model
A simulated phenotype data set is used to demonstrate model fitting. sample.col identifies the phenotype column that is matched to sample IDs in the GDS file; samples without complete model variables or without matching genotypes are excluded.
set.seed(1000)
sampid <- seqGetData(grm_fn, "sample.id") # sample IDs in the genotype file
pheno <- data.frame(
sample.id = sampid,
y = sample(c(0, 1), length(sampid), replace=TRUE, prob=c(0.95, 0.05)),
x1 = rnorm(length(sampid)),
x2 = rnorm(length(sampid)),
stringsAsFactors = FALSE)
head(pheno)## sample.id y x1 x2
## 1 HG00096 0 0.02329127 -0.1700062
## 2 HG00097 0 -1.38629108 -1.0058530
## 3 HG00099 0 0.56339700 -0.6989284
## 4 HG00100 0 -0.25703246 -1.5884806
## 5 HG00101 0 0.02926429 0.7645377
## 6 HG00102 0 -0.44480750 1.9964501grm_fn## [1] "grm_geno.gds"Null-model fitting is the computationally intensive step. The call is shown below but is not run while building this vignette:
# null model fitting using GRM from grm_fn
glmm <- seqFitNullGLMM_SPA(y ~ x1 + x2, pheno, grm_fn, trait.type="binary", sample.col="sample.id")For a fast and reproducible vignette build, the package includes a model precomputed with the same workflow:
Check the sample count, convergence status and estimated variance components before running association tests:
length(glmm$sample.id)## [1] 1092glmm$converged## [1] TRUEglmm$tau## Sigma_E Sigma_G
## 1 0Note what the components of a binary null model mean:
| Component | Meaning for trait.type="binary" |
|---|---|
coefficients |
fixed effects on the log-odds scale, including the intercept |
tau |
Sigma_E is fixed at 1, and Sigma_G is the variance of the genetic random effect on the log-odds scale |
fitted.values |
\(\mu_i = \Pr(y_i=1)\), the fitted probability of the outcome |
residuals |
\(y_i - \mu_i\), the raw residuals |
var.ratio |
variance ratio(s) used to approximate the score-statistic variance in association tests |
3.3 P-value calculations
# genetic variants stored in the file geno_fn
geno_fn## [1] "/Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/library/SeqArray/extdata/1KG_phase1_release_v3_chr22.gds"# use parallel > 1 in production to distribute variants across processes
assoc <- seqAssocGLMM_SPA(geno_fn, glmm, mac=10, parallel=1)## SAIGE association analysis:
## 2026-09-24 20:16:30
## open '1KG_phase1_release_v3_chr22.gds'
## trait type: binary
## # of samples: 1,092
## # of variants: 19,773
## MAF threshold: no
## MAC threshold: >= 10
## missing proportion threshold: <= 0.05
## variance ratio for approximation: 1
## genetic model: additive## Warning: For more accurate model building, the null model should be built using
## SAIGEgds>=v1.9.1!## # of processes: 1
##
[..................................................] 0%, ETC: --- (1/1)
[==================================================] 100%, used 1s (1/1)
[==================================================] 100%, complete, 1s
## # of variants after filtering by MAF, MAC and missing thresholds: 9,371
## P-value:
## [0,5e-10] (5e-10,5e-08] (5e-08,5e-06] (5e-06,0.0005] (0.0005,1]
## 0 0 0 3 9368
## 2026-09-24 20:16:31
## Done.head(assoc)## id chr pos rs.id ref alt AF.alt mac num
## 1 1 22 16051497 rs141578542 A G 0.30494505 666 1092
## 2 2 22 16059752 rs139717388 G A 0.05677656 124 1092
## 3 5 22 16060995 rs2843244 G A 0.06135531 134 1092
## 4 8 22 16166919 ATATTTTCTGCACATATT A 0.01190476 26 1092
## 5 9 22 16173887 GT G 0.03067766 67 1092
## 6 10 22 16205515 rs144309057 G A 0.01098901 24 1092
## beta SE pval method p.norm converged
## 1 0.2214607 0.2194430 0.3128813 Normal 0.3128813 TRUE
## 2 0.4227084 0.4240575 0.3188526 Normal 0.3188526 TRUE
## 3 0.3044670 0.4100805 0.4578107 Normal 0.4578107 TRUE
## 4 -1.0851502 0.8804521 0.2177653 Normal 0.2177653 TRUE
## 5 -0.8195449 0.5587294 0.1424302 Normal 0.1424302 TRUE
## 6 -1.0826898 0.9340256 0.2463890 Normal 0.2463890 TRUE# basic result checks
c(n.variant=nrow(assoc), n.finite.p=sum(is.finite(assoc$pval)))## n.variant n.finite.p
## 9371 9371# filtering based on p-value
assoc[assoc$pval < 5e-4, ]## id chr pos rs.id ref alt AF.alt mac num beta
## 8422 17856 22 48507315 rs57751251 T C 0.1666667 364 1092 0.971123
## 8423 17857 22 48509092 rs7292083 C T 0.1625458 355 1092 1.008351
## 8630 18282 22 49060987 rs4925399 A G 0.5929487 889 1092 0.708757
## SE pval method p.norm converged
## 8422 0.2555021 1.442051e-04 SPA 5.999175e-05 TRUE
## 8423 0.2585365 9.610252e-05 SPA 3.358045e-05 TRUE
## 8630 0.1956659 2.920160e-04 SPA 2.363433e-04 TRUEThe main result columns report the variant location and allele frequency, effect estimate (beta), standard error (SE) and p-value (pval). The method and converged columns indicate which p-value approximation was used and whether its numerical calculation converged.
The output can also be saved directly to an R object file or a GDS file:
# save to 'assoc.gds'
seqAssocGLMM_SPA(geno_fn, glmm, mac=10, parallel=1, res.savefn="assoc.gds")## SAIGE association analysis:
## 2026-09-24 20:16:32
## open '1KG_phase1_release_v3_chr22.gds'
## trait type: binary
## # of samples: 1,092
## # of variants: 19,773
## MAF threshold: no
## MAC threshold: >= 10
## missing proportion threshold: <= 0.05
## variance ratio for approximation: 1
## genetic model: additive## Warning: For more accurate model building, the null model should be built using
## SAIGEgds>=v1.9.1!## # of processes: 1
## Save to 'assoc.gds' ...
##
[..................................................] 0%, ETC: --- (1/1)
[==================================================] 100%, used 1s (1/1)
[==================================================] 100%, complete, 1s
## # of variants after filtering by MAF, MAC and missing thresholds: 9,371
## P-value:
## [0,5e-10] (5e-10,5e-08] (5e-08,5e-06] (5e-06,0.0005] (0.0005,1]
## 0 0 0 3 9368
## 2026-09-24 20:16:33
## Done.Open the output GDS file using the functions in the gdsfmt package:
# open the GDS file
(f <- openfn.gds("assoc.gds"))## File: /private/tmp/Rtmpxt1SL1/Rbuild26dc58fadfe6/SAIGEgds/vignettes/assoc.gds (506.6K)
## + [ ] *
## |--+ sample.id { Str8 1092 ZIP(26.0%), 2.2K }
## |--+ id { Int32 9371 ZIP(35.1%), 12.8K }
## |--+ chr { Str8 9371 ZIP(0.18%), 52B }
## |--+ pos { Int32 9371 ZIP(83.8%), 30.7K }
## |--+ rs.id { Str8 9371 ZIP(40.7%), 38.8K }
## |--+ ref { Str8 9371 ZIP(28.0%), 105.4K }
## |--+ alt { Str8 9371 ZIP(22.4%), 4.2K }
## |--+ AF.alt { Float64 9371 ZIP(27.4%), 20.1K }
## |--+ mac { Int32 9371 ZIP(41.6%), 15.2K }
## |--+ num { Int32 9371 ZIP(0.17%), 64B }
## |--+ beta { Float64 9371 ZIP(94.3%), 69.1K }
## |--+ SE { Float64 9371 ZIP(92.8%), 67.9K }
## |--+ pval { Float64 9371 ZIP(92.9%), 68.0K }
## |--+ method { Int32,factor 9371 ZIP(1.62%), 606B } *
## |--+ p.norm { Float64 9371 ZIP(92.9%), 68.0K }
## \--+ converged { Int32,logical 9371 ZIP(0.17%), 62B } *# get p-values
pval <- read.gdsn(index.gdsn(f, "pval"))
summary(pval)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000961 0.2862278 0.4639086 0.4980417 0.7256674 0.9999449closefn.gds(f)Load association results using the function seqSAIGE_LoadPval() in SAIGEgds:
res <- seqSAIGE_LoadPval("assoc.gds")## Loading 'assoc.gds' ...head(res)## id chr pos rs.id ref alt AF.alt mac num
## 1 1 22 16051497 rs141578542 A G 0.30494505 666 1092
## 2 2 22 16059752 rs139717388 G A 0.05677656 124 1092
## 3 5 22 16060995 rs2843244 G A 0.06135531 134 1092
## 4 8 22 16166919 ATATTTTCTGCACATATT A 0.01190476 26 1092
## 5 9 22 16173887 GT G 0.03067766 67 1092
## 6 10 22 16205515 rs144309057 G A 0.01098901 24 1092
## beta SE pval method p.norm converged
## 1 0.2214607 0.2194430 0.3128813 Normal 0.3128813 TRUE
## 2 0.4227084 0.4240575 0.3188526 Normal 0.3188526 TRUE
## 3 0.3044670 0.4100805 0.4578107 Normal 0.4578107 TRUE
## 4 -1.0851502 0.8804521 0.2177653 Normal 0.2177653 TRUE
## 5 -0.8195449 0.5587294 0.1424302 Normal 0.1424302 TRUE
## 6 -1.0826898 0.9340256 0.2463890 Normal 0.2463890 TRUE3.4 Manhattan and QQ plots for p-values
library(ggmanh)## Loading required package: ggplot2g <- manhattan_plot(assoc, pval.colname="pval", chr.colname="chr",
pos.colname="pos", x.label="Chromosome 22")
g
# QQ plot
qqunif(assoc$pval)## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggmanh package.
## Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
4 GPU Acceleration (Optional)
SAIGEgds supports optional OpenCL-based GPU acceleration for the GRM cross-product computation during null-model fitting. This can significantly speed up seqFitNullGLMM_SPA() for large sample sizes with dense genotypes.
Requirements:
- An OpenCL-compatible GPU (NVIDIA, AMD, Intel, or Apple)
- OpenCL runtime library installed on the system (no compile-time dependency)
To enable GPU acceleration, set use.gpu=TRUE:
glmm <- seqFitNullGLMM_SPA(y ~ x1 + x2, pheno, grm_fn,
trait.type="binary", sample.col="sample.id", use.gpu=TRUE)If no compatible GPU or OpenCL runtime is found, the function automatically falls back to CPU computation. When GPU initialization succeeds, use.gpu=TRUE automatically selects dense 2-bit packed genotype storage; users do not need to set geno.sparse=FALSE separately. The default sparse genotype path remains CPU-only and may be preferable when genotypes are sparse.
GPU Memory Usage:
The GPU memory required depends on the operation and dataset size. SAIGEgds checks that the total allocation does not exceed 80% of device global memory, and falls back to CPU if it does.
- Cross-product (null-model fitting): \(\lceil N/4 \rceil M\) bytes (packed genotypes) + \(16M\) bytes (standardized genotype lookup) + \(8N + 4M\) bytes (input, output and dot-product vectors), where \(N\) is the number of samples and \(M\) is the number of variants. For example, \(N=400\mathrm{K}\) and \(M=100\mathrm{K}\) requires approximately 10 GB.
- Dense GRM GEMV: \(4N^2\) bytes (GRM matrix in single precision) + \(8N\) bytes (input/output vectors). For example, \(N=5\mathrm{K}\) requires approximately 100 MB.
- GRM construction (sparse/dense): \(\lceil M/4 \rceil N\) bytes (packed genotypes) + \(16M\) bytes (standardized genotype lookup) + \(4B^2\) bytes (an implementation-defined block output buffer of width \(B\)). For \(N=400\mathrm{K}\) and \(M=2\mathrm{K}\) (the default
nsnp.sub.randomsubset), the packed genotypes require approximately 200 MB.
GPU compute buffers use single-precision (32-bit) floating point. Results are converted back to double precision on the host.
5 Session Information
sessionInfo()## R version 4.6.1 Patched (2026-06-24 r90190)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.7.8
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggmanh_1.17.0 ggplot2_4.0.3 SNPRelate_1.47.0 SAIGEgds_2.13.3
## [5] Rcpp_1.1.2 SeqArray_1.53.3 gdsfmt_1.49.8 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] tidyr_1.3.2 sass_0.4.10 generics_0.1.4
## [4] lattice_0.23-1 digest_0.6.39 magrittr_2.0.5
## [7] evaluate_1.0.5 grid_4.6.1 RColorBrewer_1.1-3
## [10] bookdown_0.48 CompQuadForm_1.4.4 fastmap_1.2.0
## [13] jsonlite_2.0.0 Matrix_1.7-6 tinytex_0.61
## [16] BiocManager_1.30.27 purrr_1.2.2 scales_1.4.0
## [19] RhpcBLASctl_0.23-42 Biostrings_2.81.9 jquerylib_0.1.4
## [22] cli_3.6.6 rlang_1.3.0 crayon_1.5.3
## [25] XVector_0.53.0 withr_3.0.3 cachem_1.1.0
## [28] yaml_2.3.12 otel_0.2.0 tools_4.6.1
## [31] parallel_4.6.1 dplyr_1.2.1 BiocGenerics_0.59.12
## [34] vctrs_0.7.3 R6_2.6.1 magick_2.9.1
## [37] stats4_4.6.1 lifecycle_1.0.5 Seqinfo_1.3.2
## [40] S4Vectors_0.51.10 IRanges_2.47.5 pkgconfig_2.0.3
## [43] pillar_1.11.1 RcppParallel_6.2.1 bslib_0.12.0
## [46] gtable_0.3.6 glue_1.8.1 tidyselect_1.2.1
## [49] tibble_3.3.1 xfun_0.61 GenomicRanges_1.65.4
## [52] knitr_1.52 dichromat_2.0-1 farver_2.1.2
## [55] htmltools_0.5.9 labeling_0.4.3 rmarkdown_2.32
## [58] compiler_4.6.1 S7_0.2.26 References
Appendix
- Zheng X, Davis JW. SAIGEgds – an efficient statistical tool for large-scale PheWAS with mixed models. Bioinformatics (2021). Mar;37(5):728-730. DOI: 10.1093/bioinformatics/btaa731.
- Zhou W, Nielsen JB, Fritsche LG, Dey R, Gabrielsen ME, Wolford BN, LeFaive J, VandeHaar P, Gagliano SA, Gifford A, Bastarache LA, Wei WQ, Denny JC, Lin M, Hveem K, Kang HM, Abecasis GR, Willer CJ, Lee S. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nat Genet (2018). Sep;50(9):1335-1341. DOI: 10.1038/s41588-018-0184-y
- Zheng X, Gogarten S, Lawrence M, Stilp A, Conomos M, Weir BS, Laurie C, Levine D. SeqArray – A storage-efficient high-performance data format for WGS variant calls. Bioinformatics (2017). DOI: 10.1093/bioinformatics/btx145.
A See also
Time-to-event outcomes are covered in Time-to-event (survival) analysis in SAIGEgds (vignette("SAIGEgds_surv", package="SAIGEgds")). Rare-variant set-based tests using full and sparse GRMs are covered in Set-based Rare-Variant Association Analysis in SAIGEgds (vignette("SAIGEgds_set-aggr", package="SAIGEgds")).
SeqArray: Data Management of Large-scale Whole-genome Sequence Variant Calls
SNPRelate: Parallel Computing Toolset for Relatedness and Principal Component Analysis of SNP Data