GPTomics/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

| Method | Citation | Mechanism | When | |--------|----------|-----------|------| | Burden / collapsing (CMC, weighted-sum) | Li & Leal 2008; Madsen & Browning 2009 | Collapse variants into one score, test its single coefficient; assumes one effect direction | Strong prior that variants act the same way (e.g. LoF in a gene) | | SKAT | Wu 2011 | Variance-component score test on summed squared weighted single-variant scores; directions do not cancel | Mixed directions, or many neutral variants diluting the set | | SKAT-O | Lee 2012 | Optimal linear combination of burden and SKAT over a rho grid in [0,1]; data choose rho | Unknown architecture - the safe default | | ACAT-V / ACAT-O | Liu 2019 | Cauchy combination of p-values, calibrated under arbitrary dependence, no permutation/GRM; the smallest p dominates (one artifact can drive it) | Sparse-causal sets, fast omnibus, combining masks/tests | | STAAR / STAAR-O | Li 2020 | Variance-component test weighting variants by multiple functional annotations (annotation PCs) | WGS regulatory regions where annotations carry the signal | | SAIGE-GENE+ | Zhou 2022 | LMM null + SPA + variance ratio; multiple MAF cutoffs and annotations in one set test | Biobank binary traits, case/control imbalance, relatedness | | regenie --vc-tests | Mbatchou 2021 | Whole-genome ridge null (step 1), then masked burden/SKAT/SKAT-O/ACAT in step 2 with Firth/SPA | Biobank pipelines wanting single-variant and gene tests together |

Decision Tree by Scenario

| Scenario | Use | Why | |----------|-----|-----| | Strong prior all variants act one direction (LoF mask) | Burden / collapsing | Most powerful under a true single direction | | Risk and protective variants expected in the same set | SKAT | Squared scores, directions do not cancel | | Architecture unknown | SKAT-O | Optimizes rho between burden and SKAT | | Sparse causal set, or one omnibus across masks | ACAT-V / ACAT-O | Dependence-robust Cauchy combiner, no permutation | | WGS noncoding where annotations carry the signal | STAAR-O | Multiple functional-annotation weights in one test | | Biobank, imbalanced binary trait, relatedness | SAIGE-GENE+ | SPA + LMM null keeps the tail calibrated at low MAC | | One pipeline for single-variant + gene tests at biobank scale | regenie --vc-tests | Shared step-1 null, Firth/SPA in step 2 | | Small cohort, full control of mask and weights | SKAT R package | Direct, scriptable, SSD files for many sets |

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).

SKAT underpowered when truth is unidirectional

Trigger: a clean LoF mask where every variant raises risk. Mechanism: the variance-component test spends power on a 2-sided alternative it does not need. Symptom: burden hits, SKAT does not. Fix: SKAT-O (lets rho->1) or a burden test for that mask.

Wrong or default mask

Trigger: running one default MAF cutoff or an unfiltered annotation set. Mechanism: the mask defines the hypothesis; a too-loose MAF or LoF+benign-missense mask dilutes signal with noise. Symptom: a true gene disappears under one mask, appears under another. Fix: test a small grid of masks (LoF, LoF+missense; MAF 0.001, 0.01) and combine with ACAT-O, accounting for the masks in the burden.

Uncalibrated tail under case/control imbalance

Trigger: a naive set test on an imbalanced binary trait at low MAC. Mechanism: the score statistic's normal/chi-square null is wrong in the tail. Symptom: anti-conservative gene p-values, inflated QQ for rare masks. Fix: SAIGE-GENE+ (SPA) or regenie --firth/--spa; never a plain score test here.

Residual structure/relatedness in the null

Trigger: aggregating in a structured or related sample with a fixed-effect-only null. Mechanism: relatedness is a covariance structure PCs cannot remove. Symptom: genome-wide gene inflation that PCs do not fix. Fix: an LMM null (SAIGE-GENE+, regenie step-1 LOCO predictor) before the set test.

Imputed/low-quality variants in the mask

Trigger: building masks from imputed or low-callrate genotypes. Mechanism: miscalled rare variants add spurious carriers. Symptom: unreplicable gene hits driven by a few low-quality sites. Fix: filter by INFO/R2 (>=0.8 for rare) and genotype quality before masking; prefer sequenced calls for rare-variant sets.

Quantitative Thresholds

| Quantity | Typical value | Rationale | |----------|---------------|-----------| | "Rare" MAF cutoff for aggregation | MAF < 0.01 (< 0.001 for LoF-only) | below this single-variant power collapses; aggregate instead | | Mask MAF tiers | 0.0001 / 0.001 / 0.01 | nested AAF bins capture ultra-rare and rare jointly (regenie --aaf-bins, SAIGE --maxMAF_in_groupTest) | | Exome-wide gene significance | ~2.5e-6 | Bonferroni 0.05 over ~20,000 genes; tighten further for multiple masks per gene | | Per-mask multiple testing | divide by (genes x masks) or combine masks via ACAT-O | each mask is a separate test unless an omnibus absorbs them | | Ultra-rare collapsing (MAC) | collapse MAC < ~10 into one pseudo-variant | regenie --vc-MACthr default 10; SAIGE-GENE+ collapses ultra-rare for calibration | | Imputed-variant quality for masks | INFO/R2 >= 0.8 (rare) | rare imputed dosages are noisy; lenient 0.3 cutoffs corrupt masks | | SKAT MAF weighting | Beta(MAF; 1, 25) | weights.beta=c(1,25) up-weights rarer variants (Wu 2011 default) |

Thresholds are conventions; inspect the per-gene QQ plot and verify current best practice before applying numbers blindly.

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | regenie --vc-tests skat-o unrecognized | wrong token | use skato (and acato, acatv, skato-acat); not skat-o | | Empty or tiny masks, genes silently dropped | variants in set-list absent from anno-file | run --check-burden-files first; harmonize variant IDs | | Inflated gene QQ for an imbalanced trait | plain score test, no SPA/Firth | SAIGE-GENE+ or regenie --firth/--spa | | SKAT-O not actually run in R | passing method="SKAT" or omitting it | method="SKATO" or "optimal.adj"; r.corr=1 is pure burden | | Single MAF cutoff misses ultra-rare signal | one --aaf-bins/--maxMAF value | pass nested cutoffs (0.0001,0.001,0.01) in one run | | Same single-variant p reported as "gene" | testing markers, not a set | confirm a set/group file is supplied and the test is set-based | | Gene hit driven by one artifactual variant | ACAT/burden dominated by a miscalled site | QC inputs (INFO/R2, genotype quality) before aggregating |

References

1. Li B, Leal SM. Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. American Journal of Human Genetics 2008; 83(3):311-321. DOI:10.1016/j.ajhg.2008.06.024.
2. Madsen BE, Browning SR. A groupwise association test for rare mutations using a weighted sum statistic. PLoS Genetics 2009; 5(2):e1000384. DOI:10.1371/journal.pgen.1000384.
3. Wu MC, Lee S, Cai T, Li Y, Boehnke M, Lin X. Rare-variant association testing for sequencing data with the sequence kernel association test. American Journal of Human Genetics 2011; 89(1):82-93. DOI:10.1016/j.ajhg.2011.05.029.
4. Lee S, Emond MJ, Bamshad MJ, et al. Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies. American Journal of Human Genetics 2012; 91(2):224-237. DOI:10.1016/j.ajhg.2012.06.007.
5. 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. American Journal of Human Genetics 2019; 104(3):410-421. DOI:10.1016/j.ajhg.2019.01.002.
6. Li X, Li Z, Zhou H, et al. Dynamic incorporation of multiple in silico functional annotations empowers rare variant association analysis of large whole-genome sequencing studies at scale. Nature Genetics 2020; 52(9):969-983. DOI:10.1038/s41588-020-0676-4.
7. Mbatchou J, Barnard L, Backman J, et al. Computationally efficient whole-genome regression for quantitative and binary traits. Nature Genetics 2021; 53(7):1097-1103. DOI:10.1038/s41588-021-00870-7.
8. Zhou W, Bi W, Zhao Z, et al. SAIGE-GENE+ improves the efficiency and accuracy of set-based rare variant association tests. Nature Genetics 2022; 54(10):1466-1469. DOI:10.1038/s41588-022-01178-w.

Related Skills

• association-testing - single-variant GWAS (linear/logistic/LMM/SPA per marker) that this skill aggregates beyond
• plink-basics - genotype QC and format conversion before masking
• population-structure - PCs and relatedness for the null model that calibrates the set test
• variant-calling/variant-annotation - functional annotations (LoF, missense, CADD, regulatory) that define masks
• clinical-databases/variant-prioritization - clinical interpretation of variants flagged by gene tests
• causal-genomics/fine-mapping - resolving which variants in a significant gene/region carry the signal