Bio Population Genetics Rare Variant Association

Version Compatibility

Reference examples tested with: regenie 3.4+, SAIGE 1.3+, SKAT 2.2+ (R), STAAR 0.9.7+ (R).

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Version traps that change results, not just syntax: regenie --vc-tests accepts skat,skato,skato-acat,acatv,acato,acato-full (not skat-o), --aaf-bins upper bounds always add an implicit singleton mask, and --build-mask defaults to max (one carrier-status column per set) not sum. SAIGE-GENE+ --maxMAF_in_groupTest takes multiple comma-separated cutoffs in ONE run (the whole point of GENE+ over GENE). The SKAT R package selects SKAT-O with method="optimal.adj" or method="SKATO", and r.corr is the rho grid (0=SKAT, 1=burden); weights.beta=c(1,25) is the rarer-up-weighting default. The single source of truth for versions is this block, not headings.

Rare-Variant Association

"Test whether rare variants in this gene associate with my trait" -> Aggregate the rare variants in a gene or region into one set-based statistic under an explicit mask, because no single rare variant has enough carriers to test alone.

• CLI: regenie --step 2 --anno-file ... --set-list ... --mask-def ... --aaf-bins 0.01 --vc-tests skato,acato (biobank masks plus omnibus tests)
• CLI: step2_SPAtests.R --groupFile ... --annotation_in_groupTest lof,missense;lof --maxMAF_in_groupTest 0.0001,0.001,0.01 (SAIGE-GENE+, imbalance-robust)
• R: SKAT(Z, obj, method="SKATO", weights.beta=c(1,25)) (direct, small cohorts)

Scope: gene/region-based rare-variant aggregation (burden, SKAT, SKAT-O, ACAT-V/ACAT-O, STAAR), variant masks (functional class plus MAF cutoff), and the per-gene multiple-testing burden. Single-variant GWAS (linear/logistic/LMM/SPA per marker) routes to association-testing. The functional annotations that define masks (LoF, missense, CADD, regulatory) come from variant-calling/variant-annotation. Variant prioritization for clinical interpretation routes to clinical-databases/variant-prioritization.

The Single Most Important Insight -- a gene-based test is a bet about the direction-of-effect architecture, and the mask IS the hypothesis

1. Single-variant GWAS is underpowered for rare variants because a handful of carriers gives a tiny non-centrality, so signal must be aggregated across a gene or region - and HOW it is aggregated encodes a belief about the unobserved effect architecture.
2. A burden test collapses the set into one direction and is the most powerful test WHEN that holds, but mixing risk and protective variants makes their contributions cancel to a null (false negative); SKAT sums squared scores so directions cannot cancel but is weaker when the truth is unidirectional; SKAT-O optimizes a rho grid between the two and is the default when the architecture is unknown.
3. The MASK is the hypothesis, not a preprocessing detail: which variants enter (LoF-only vs LoF+missense, MAF<0.01 vs <0.001, annotation weights) defines what is being tested, and a different mask is a different question with a different answer - so report the mask, not just the p-value.
4. The aggregate is only as calibrated as the null model: case/control imbalance and low MAC make naive set tests anti-conservative (SAIGE-GENE+ uses SPA, regenie uses Firth/SPA), population structure and relatedness still need an LMM null, and imputed/low-quality variants silently corrupt the mask unless filtered by INFO/R2 and genotype quality first.

Tool Taxonomy

Decision Tree by Scenario

Build the Mask and Run Aggregate Tests with regenie

Goal: test each gene under one or more masks (functional class x MAF cutoff) using burden plus variance-component tests in a biobank-scale pipeline.

Approach: reuse the step-1 whole-genome ridge null, then in step 2 define annotations, gene sets, and mask rules and request the omnibus tests, letting Firth handle imbalanced binary traits.

# Step 1 builds the LOCO whole-genome predictor (the null) once, shared with single-variant GWAS.
regenie --step 1 --bed geno_array --phenoFile pheno.txt --covarFile covar.txt \
    --bsize 1000 --lowmem --out fit_null

# Step 2: --anno-file maps variant -> gene -> annotation; --set-list lists each gene's variants;
# --mask-def names which annotation categories form each mask. --aaf-bins sets the MAF ceilings
# (a singleton mask is always added). --vc-tests requests SKAT-O and the ACAT omnibus alongside
# burden. --firth keeps the imbalanced binary-trait tail calibrated; --build-mask max is the default.
regenie --step 2 --bed geno_wes --phenoFile pheno.txt --covarFile covar.txt \
    --pred fit_null_pred.list --anno-file annot.txt --set-list sets.txt --mask-def masks.txt \
    --aaf-bins 0.001,0.01 --vc-tests skato,acato --build-mask max \
    --bt --firth --approx --pThresh 0.05 --out gene_tests

Mask-building file formats (one entry per line, space/tab separated):

• annot.txt: VARIANT_ID GENE ANNOTATION (e.g. 1:55039839:T:C PCSK9 LoF); variants with no entry fall in NULL.
• sets.txt: GENE CHR POS VARIANT_ID,VARIANT_ID,... (the gene plus its comma-separated variant list).
• masks.txt: MASK_NAME ANNOTATION,ANNOTATION (e.g. Mask_LoF LoF and Mask_LoF_mis LoF,missense).

Run regenie --step 2 ... --check-burden-files --ignore-pred first to catch variants in the set-list that are absent from the annotation file (a silent source of empty or wrong masks). --ignore-pred is required here because this validation runs before the step-1 predictor exists.

Imbalance-Robust Set Tests with SAIGE-GENE+

Goal: test genes for an imbalanced binary trait in a related sample, scanning several MAF cutoffs and annotation groups in one pass.

Approach: fit the SPA-LMM null once (step 1, with a variance ratio), then run the set test passing multiple annotations and multiple max-MAF thresholds so GENE+ combines them.

# Step 2 set test. --annotation_in_groupTest gives the masks (semicolon-separated groups, each a
# comma-separated annotation list). --maxMAF_in_groupTest passes several MAF cutoffs in ONE run -
# this multi-cutoff combination is exactly what GENE+ adds over the original SAIGE-GENE.
step2_SPAtests.R --bgenFile geno_wes.bgen --groupFile groups.txt \
    --GMMATmodelFile null.rda --varianceRatioFile null.varianceRatio.txt \
    --annotation_in_groupTest "lof;lof,missense;lof,missense,synonymous" \
    --maxMAF_in_groupTest 0.0001,0.001,0.01 --is_output_moreDetails TRUE \
    --SAIGEOutputFile gene_tests.txt

The groups.txt file gives, per gene, a line of variant IDs and a matching line of their annotations (and optionally a weight line); the annotation labels there must match --annotation_in_groupTest.

Direct SKAT-O in R for a Small Cohort

Goal: run burden, SKAT, and SKAT-O on a gene's rare-variant genotype matrix with explicit MAF weighting, for a sample small enough to hold in memory.

Approach: fit the null model once on covariates, then call SKAT per gene with the rho grid; for many genes use SSD files keyed by a SetID rather than passing matrices.

library(SKAT)

# Null model on covariates only (out_type='D' binary, 'C' continuous). Refit once, reuse per gene.
obj <- SKAT_Null_Model(phenotype ~ age + sex + PC1 + PC2, out_type = 'D', data = covar_df)

# Z is the n x m genotype matrix (0/1/2) for the m rare variants in one gene.
# weights.beta=c(1,25) is the Beta(MAF;1,25) up-weighting of rarer variants (the SKAT default; the
# Madsen-Browning weight is the gentler Beta(0.5,0.5)). method='SKATO' searches the rho grid (rho=0
# SKAT, rho=1 burden); 'burden' or default SKAT recover the endpoints. r.corr passes an explicit rho grid.
skato <- SKAT(Z, obj, method = 'SKATO', weights.beta = c(1, 25))
burden <- SKAT(Z, obj, r.corr = 1, weights.beta = c(1, 25))
skat <- SKAT(Z, obj, weights.beta = c(1, 25))
c(skato = skato$p.value, burden = burden$p.value, skat = skat$p.value)

For genome-wide gene scans, build an SSD file with Generate_SSD_SetID(bed, bim, fam, SetID, SSD, Info), Open_SSD(), then SKAT.SSD.All(SSD.INFO, obj) to test every set without holding all matrices in memory.

Per-Method Failure Modes

Burden test cancels under mixed directions

Trigger: a mask mixing risk and protective (or many null) variants. Mechanism: the collapsed score sums signed contributions that offset. Symptom: near-null p for a gene that SKAT flags strongly. Fix: use SKAT or SKAT-O; reserve pure burden for a mask with a real single-direction prior (LoF-only).