Set-based Rare-Variant Association Analysis in SAIGEgds

Introduction

Set-based tests aggregate rare variants within a gene or genomic region. They can gain power over single-variant tests when several variants in the set affect the same phenotype, but low minor allele counts make the tests particularly sensitive to sample relatedness. SAIGEgds implements Burden, SKAT, ACAT-V and ACAT-O tests in the SAIGE/SAIGE-GENE framework [1,2].

For set-based analysis, two complementary genetic relationship matrices (GRMs) are used in the same null model:

  • the full GRM, represented by a set of common, preferably LD-pruned markers, models the genome-wide polygenic effect and population/sample structure during null-model fitting;
  • the sparse GRM retains the diagonal and related sample pairs above a chosen relatedness threshold. It supplies the covariance projection needed by SKAT and ACAT-O and improves the calibration of rare-variant statistics in samples that include relatives.

The sparse GRM supplements rather than replaces the full GRM. Its main benefit is more accurate modeling of rare-variant score covariance; it does not guarantee a power increase for every phenotype or variant set.

This vignette demonstrates a binary-trait workflow. Quantitative traits use the same functions with trait.type="quantitative". For survival outcomes, only Burden and ACAT-V are currently available; see the companion survival vignette.

Input data

Three inputs are kept separate in production analyses:

  1. phenotype and covariate data with one row per sample;
  2. common, LD-pruned markers for the full and sparse GRMs;
  3. rare variants to be grouped and tested, usually from exome or whole-genome sequencing data.

The sample IDs must be shared across these inputs. Their original orders do not need to match because SAIGEgds matches samples by ID.

library(Matrix)
library(SeqArray)
library(SAIGEgds)

# Common markers used by the GRMs: 1,000 samples and 10,000 SNPs
grm.fn <- system.file("extdata", "grm1k_10k_snp.gds", package="SAIGEgds")

# For this compact example, variants to group and test use the same file.
# Production analyses should use a separate exome or whole-genome GDS file.
assoc.fn <- grm.fn

# Binary phenotype and covariates
pheno.fn <- system.file("extdata", "pheno.txt.gz", package="SAIGEgds")
pheno <- read.table(pheno.fn, header=TRUE, as.is=TRUE)
head(pheno)
##   sample.id y     yy      x1 x2
## 1        s1 0 4.5542  1.5118  1
## 2        s2 0 3.7941  0.3898  1
## 3        s3 0 5.0411 -0.6212  1
## 4        s4 0 5.6394 -2.2147  1
## 5        s5 0 4.2134  1.1249  1
## 6        s6 0 4.6145 -0.0449  1
table(pheno$y)
## 
##   0   1 
## 902  98

Constructing the two GRMs

Full GRM

The full GRM does not need to be materialized as an \(N \times N\) matrix. Passing grm.fn to seqFitNullGLMM_SPA() loads the standardized marker genotypes and computes full-GRM matrix-vector products internally. In a real analysis, the GRM marker file should contain high-quality common variants pruned for linkage disequilibrium across the genome.

Sparse GRM

seqFitSparseGRM() first uses a random marker subset to identify potentially related pairs and then recalculates those entries using the full candidate marker set. With the default rel.cutoff=0.125, entries below the threshold are omitted apart from the diagonal.

grm.gds <- seqOpen(grm.fn)
sp.grm <- seqFitSparseGRM(grm.gds, sample.id=pheno$sample.id,
    nsnp.sub.random=2000, rel.cutoff=0.125, num.thread=2)
seqClose(grm.gds)

# Save once and reuse for phenotypes measured in the same cohort
saveRDS(sp.grm, "sparse_grm.rds")

The package includes a sparse GRM precomputed by this workflow so the vignette does not rebuild it each time.

sp.grm.fn <- system.file("extdata", "grm1k_10k_sp_grm.rds",
    package="SAIGEgds")
sp.grm <- readRDS(sp.grm.fn)

dim(sp.grm)
## [1] 1000 1000
class(sp.grm)
## [1] "dsCMatrix"
## attr(,"package")
## [1] "Matrix"
nnzero(sp.grm)
## [1] 6642
nnzero(sp.grm) / prod(dim(sp.grm))
## [1] 0.006642
stopifnot(setequal(colnames(sp.grm), pheno$sample.id))

For a one-step analysis, grm.mat=TRUE in seqFitNullGLMM_SPA() calls seqFitSparseGRM() internally. Constructing the sparse GRM explicitly is usually preferable when many phenotypes share the same samples because the matrix can be saved and reused.

Fitting the null model

Supply both GRMs in one call: gdsfile=grm.fn provides the full GRM and grm.mat=sp.grm provides the sparse GRM. use.cateMAC=TRUE estimates variance ratios by minor allele count (MAC) category instead of applying one global ratio; this is useful when the tested variants span ultra-rare through low-frequency categories.

glmm <- seqFitNullGLMM_SPA(y ~ x1 + x2, pheno,
    gdsfile=grm.fn, grm.mat=sp.grm,
    trait.type="binary", sample.col="sample.id",
    use.cateMAC=TRUE, num.thread=2, verbose=FALSE)

Check convergence and confirm that the sparse-GRM projection components were stored. SKAT and ACAT-O require these components.

glmm$converged
## [1] TRUE
glmm$tau
##   Sigma_E   Sigma_G 
## 1.0000000 0.3368194
head(glmm$var.ratio)
##         id    maf mac       var1     var2_u     ratio
## 1 simu1_g1 0.0010   2 0.04374719 0.04375442 0.9998349
## 2 simu1_g2 0.0020   4 0.08146448 0.08146984 0.9999342
## 3 simu1_g3 0.0010   2 0.09156725 0.09158307 0.9998272
## 4 simu1_g4 0.0015   3 0.04801201 0.04802110 0.9998107
## 5 simu1_g5 0.0020   4 0.08466725 0.08464713 1.0002376
## 6 simu1_g6 0.0025   5 0.08584211 0.08581817 1.0002790
c(Sigma_inv = !is.null(glmm$Sigma_inv),
    chol_inv_X_Sigma = !is.null(glmm$chol_inv_X_Sigma))
##        Sigma_inv chol_inv_X_Sigma 
##             TRUE             TRUE

The equivalent one-step call below derives the sparse GRM from the full-GRM marker file during fitting:

glmm <- seqFitNullGLMM_SPA(y ~ x1 + x2, pheno,
    gdsfile=grm.fn, grm.mat=TRUE,
    trait.type="binary", sample.col="sample.id",
    use.cateMAC=TRUE, nsnp.sub.random=2000, rel.cutoff=0.125)

Defining variant sets

The aggregate-test functions take a SeqUnitListClass object. Here, sliding windows over the bundled genotype file provide a compact reproducible example. In a gene-based analysis, use a separate exome or whole-genome association GDS file and create the same type of object from gene boundaries and functional annotations, with one unit per gene or annotation mask.

assoc.gds <- seqOpen(assoc.fn)
units <- seqUnitSlidingWindows(assoc.gds, win.size=5000, win.shift=5000)
## Chromosome 1, # of units: 2
## Chromosome 2, # of units: 2
## # of units in total: 4
units
## List of 2
##  $ desp :'data.frame':   4 obs. of  3 variables:
##   ..$ chr  : chr [1:4] "1" "1" "2" "2"
##   ..$ start: int [1:4] 0 5000 5000 10000
##   ..$ end  : int [1:4] 4999 9999 9999 14999
##  $ index:List of 4
##   ..$ chr1: int [1:4999] 1 2 3 4 5 6 7 8 9 10 ...
##   ..$ chr1: int [1:4989] 5000 5001 5002 5003 5004 5005 5006 5007 5008 5009 ...
##   ..$ chr2: int [1:11] 9989 9990 9991 9992 9993 9994 9995 9996 9997 9998 ...
##   ..$ chr2: int 10000
##  - attr(*, "class")= chr "SeqUnitListClass"

Variants are filtered within each set by the maxMAF and missing arguments of the association functions. Multiple MAF thresholds can be tested in one call. The default beta weights include Beta(1,1) and Beta(1,25), the latter upweighting rarer variants. Variants with MAC <= collapse.mac are collapsed before SKAT, ACAT-V and ACAT-O calculations.

Set-based association tests

Burden

The burden test collapses variants into a single weighted genotype score. It is most powerful when most causal effects in a set point in the same direction.

burden <- seqAssocGLMM_Burden(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), parallel=1, verbose=FALSE)
burden
##    chr start  end maxMAF numvar macmin macmed macmax summac weight        beta
## 1    1     0 4999  0.010    203      5     16     20   3126  (1,1)  0.01457060
## 2    1     0 4999  0.010    203      5     16     20   3126 (1,25)  0.01002053
## 3    1     0 4999  0.005     28      5      9     10    243  (1,1) -0.50455566
## 4    1     0 4999  0.005     28      5      9     10    243 (1,25) -0.50588254
## 5    1     0 4999  0.010      4    NaN    NaN    NaN    NaN Cauchy         NaN
## 6    1  5000 9999  0.010    209      4     17     20   3373  (1,1) -0.06743052
## 7    1  5000 9999  0.010    209      4     17     20   3373 (1,25) -0.06774527
## 8    1  5000 9999  0.005     15      4      9     10    129  (1,1) -0.13141090
## 9    1  5000 9999  0.005     15      4      9     10    129 (1,25) -0.14107367
## 10   1  5000 9999  0.010      4    NaN    NaN    NaN    NaN Cauchy         NaN
##            SE       pval method     p.norm converged
## 1  0.06206550 0.81439375 Normal 0.81439375      TRUE
## 2  0.06207190 0.87175138 Normal 0.87175138      TRUE
## 3  0.23384730 0.03095671    SPA 0.03121860      TRUE
## 4  0.23377520 0.03046664    SPA 0.03072198      TRUE
## 5         NaN 0.07438731   <NA>        NaN      TRUE
## 6  0.06141192 0.27220285 Normal 0.27220285      TRUE
## 7  0.06153948 0.27096525 Normal 0.27096525      TRUE
## 8  0.32388370 0.68493743 Normal 0.68493743      TRUE
## 9  0.32350916 0.66278363 Normal 0.66278363      TRUE
## 10        NaN 0.45823687   <NA>        NaN      TRUE

SKAT

SKAT is a variance-component test that allows effect directions and magnitudes to vary. The sparse-GRM projection in the null model is required to construct the within-set covariance matrix while accounting for related samples.

skat <- seqAssocGLMM_SKAT(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), collapse.mac=10, parallel=1, verbose=FALSE)
skat
##    chr start  end maxMAF numvar macmin macmed macmax summac weight n_collapse
## 1    1     0 4999  0.010    203      5     16     20   3126  (1,1)         28
## 2    1     0 4999  0.010    203      5     16     20   3126 (1,25)         28
## 3    1     0 4999  0.005     28      5      9     10    243  (1,1)         28
## 4    1     0 4999  0.005     28      5      9     10    243 (1,25)         28
## 5    1     0 4999  0.010      4    NaN    NaN    NaN    NaN Cauchy         NA
## 6    1  5000 9999  0.010    209      4     17     20   3373  (1,1)         15
## 7    1  5000 9999  0.010    209      4     17     20   3373 (1,25)         15
## 8    1  5000 9999  0.005     15      4      9     10    129  (1,1)         15
## 9    1  5000 9999  0.005     15      4      9     10    129 (1,25)         15
## 10   1  5000 9999  0.010      4    NaN    NaN    NaN    NaN Cauchy         NA
##    g_ncol g_minMAC       pval
## 1     176        0 0.07002189
## 2     176        0 0.38679655
## 3       1        0 0.03663650
## 4       1        0 0.03663650
## 5      NA      NaN 0.05688580
## 6     195        0 0.97183997
## 7     195        0 0.96207920
## 8       1        0 0.80619009
## 9       1        0 0.80619009
## 10     NA      NaN 0.94398905

ACAT-V

ACAT-V combines variant-level p-values using a Cauchy combination. It can be effective when only a small fraction of variants in a set are causal.

acatv <- seqAssocGLMM_ACAT_V(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), collapse.mac=10, parallel=1, verbose=FALSE)
acatv
##    chr start  end maxMAF numvar macmin macmed macmax summac weight n_single
## 1    1     0 4999  0.010    203      5     16     20   3126  (1,1)      175
## 2    1     0 4999  0.010    203      5     16     20   3126 (1,25)      175
## 3    1     0 4999  0.005     28      5      9     10    243  (1,1)        0
## 4    1     0 4999  0.005     28      5      9     10    243 (1,25)        0
## 5    1     0 4999  0.010      4    NaN    NaN    NaN    NaN Cauchy       NA
## 6    1  5000 9999  0.010    209      4     17     20   3373  (1,1)      194
## 7    1  5000 9999  0.010    209      4     17     20   3373 (1,25)      194
## 8    1  5000 9999  0.005     15      4      9     10    129  (1,1)        0
## 9    1  5000 9999  0.005     15      4      9     10    129 (1,25)        0
## 10   1  5000 9999  0.010      4    NaN    NaN    NaN    NaN Cauchy       NA
##    n_collapse       pval
## 1          28 0.13392498
## 2          28 0.13877741
## 3          28 0.03095671
## 4          28 0.03095671
## 5          NA 0.05073845
## 6          15 0.82183036
## 7          15 0.82992407
## 8          15 0.68493743
## 9          15 0.68493743
## 10         NA 0.77214222

ACAT-O

ACAT-O combines Burden, SKAT and ACAT-V results, providing an omnibus test that is robust to different genetic architectures. Because SKAT is one of its components, ACAT-O also requires the sparse GRM in the fitted null model.

acato <- seqAssocGLMM_ACAT_O(assoc.gds, glmm, units,
    maxMAF=c(0.01, 0.005), collapse.mac=10, parallel=1, verbose=FALSE)
acato
##    chr start  end maxMAF numvar macmin macmed macmax summac weight n_collapse
## 1    1     0 4999  0.010    203      5     16     20   3126  (1,1)         28
## 2    1     0 4999  0.010    203      5     16     20   3126 (1,25)         28
## 3    1     0 4999  0.005     28      5      9     10    243  (1,1)         28
## 4    1     0 4999  0.005     28      5      9     10    243 (1,25)         28
## 5    1     0 4999  0.010      4    NaN    NaN    NaN    NaN Cauchy         NA
## 6    1  5000 9999  0.010    209      4     17     20   3373  (1,1)         15
## 7    1  5000 9999  0.010    209      4     17     20   3373 (1,25)         15
## 8    1  5000 9999  0.005     15      4      9     10    129  (1,1)         15
## 9    1  5000 9999  0.005     15      4      9     10    129 (1,25)         15
## 10   1  5000 9999  0.010      4    NaN    NaN    NaN    NaN Cauchy         NA
##          pval   p.burden     p.skat    p.acatv burden.beta  burden.se
## 1  0.16680099 0.81439375 0.07002189 0.13392498  0.01457060 0.06206550
## 2  0.48180146 0.87175138 0.38679655 0.13877741  0.01002053 0.06207190
## 3  0.03264436 0.03095671 0.03663650 0.03095671 -0.50455566 0.23384730
## 4  0.03246089 0.03046664 0.03663650 0.03095671 -0.50588254 0.23377520
## 5  0.05915300        NaN        NaN        NaN         NaN        NaN
## 6  0.92202502 0.27220285 0.97183997 0.82183036 -0.06743052 0.06141192
## 7  0.89933737 0.27096525 0.96207920 0.82992407 -0.06774527 0.06153948
## 8  0.73602012 0.68493743 0.80619009 0.68493743 -0.13141090 0.32388370
## 9  0.73042711 0.66278363 0.80619009 0.68493743 -0.14107367 0.32350916
## 10 0.86496007        NaN        NaN        NaN         NaN        NaN

Choosing a test

Comparison of set-based tests.
Test Alternative captured well Sparse GRM required
Burden many effects in the same direction no
SKAT effects vary in direction or magnitude yes
ACAT-V a small fraction of variants drive the signal no
ACAT-O unknown architecture; omnibus combination yes

Although Burden and ACAT-V can run without a sparse GRM, using the same full + sparse null model for all four tests gives a consistent relatedness adjustment. Report the tested annotation mask, MAF threshold, weight, ultra-rare collapsing threshold and test name together with each p-value.

Practical recommendations

  • Build the full and sparse GRMs from common, high-quality, LD-pruned markers, not from the rare variants being tested.
  • Reuse one sparse GRM only when the analyzed sample cohort is the same. SAIGEgds subsets and reorders a named matrix, but its row and column names must be sample IDs.
  • Inspect the sparse matrix density. A very low cutoff can make it nearly dense, increasing memory and factorization time; a very high cutoff can omit meaningful relatedness.
  • Use category-specific variance ratios (use.cateMAC=TRUE) when the analysis includes very low MAC variants, and retain the default simulated markers for MAC categories with too few observed GRM markers.
  • Define biologically motivated masks before association testing. Common choices include loss-of-function alone and loss-of-function plus damaging missense variants.
  • Save large results directly with res.savefn and distribute units with the parallel argument in production runs.

Session Information

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] SAIGEgds_2.13.3  Rcpp_1.1.2       SeqArray_1.53.3  gdsfmt_1.49.8   
## [5] Matrix_1.7-6     BiocStyle_2.41.0
## 
## loaded via a namespace (and not attached):
##  [1] jsonlite_2.0.0       compiler_4.6.1       BiocManager_1.30.27 
##  [4] crayon_1.5.3         GenomicRanges_1.65.4 Biostrings_2.81.9   
##  [7] parallel_4.6.1       survey_4.5           jquerylib_0.1.4     
## [10] splines_4.6.1        IRanges_2.47.5       Seqinfo_1.3.2       
## [13] yaml_2.3.12          fastmap_1.2.0        lattice_0.23-1      
## [16] R6_2.6.1             XVector_0.53.0       generics_0.1.4      
## [19] knitr_1.52           BiocGenerics_0.59.12 maketools_1.3.2     
## [22] DBI_1.3.0            bslib_0.12.0         rlang_1.3.0         
## [25] cachem_1.1.0         xfun_0.61            sass_0.4.10         
## [28] sys_3.4.3            RcppParallel_6.2.1   cli_3.6.6           
## [31] SKAT_2.2.5           digest_0.6.39        grid_4.6.1          
## [34] SPAtest_3.1.2        CompQuadForm_1.4.4   lifecycle_1.0.5     
## [37] S4Vectors_0.51.10    RSpectra_0.16-2      evaluate_1.0.5      
## [40] mitools_2.7          buildtools_1.0.0     survival_3.8-12     
## [43] stats4_4.6.1         rmarkdown_2.32       tools_4.6.1         
## [46] htmltools_0.5.9

References

  1. Zhou W, Nielsen JB, Fritsche LG, et al. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nat Genet (2018). 50:1335-1341. DOI: 10.1038/s41588-018-0184-y.
  2. Zhou W, Zhao Z, Nielsen JB, et al. Scalable generalized linear mixed model for region-based association tests in large biobanks and cohorts. Nat Genet (2020). 52:634-639. DOI: 10.1038/s41588-020-0621-6.
  3. Liu Y, Chen S, Li Z, Morrison AC, Boerwinkle E, Lin X. ACAT: A fast and powerful p value combination method for rare-variant analysis in sequencing studies. Am J Hum Genet (2019). 104:410-421. DOI: 10.1016/j.ajhg.2019.01.002.

See also

Single-variant tests are covered in SAIGEgds Tutorial (single variant tests) (vignette("SAIGEgds", package="SAIGEgds")). Time-to-event outcomes are covered in SAIGEgds Tutorial (time-to-event outcomes) (vignette("SAIGEgds_surv", 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