Gene-Based Burden Tests for Rare Variants: From Foundational Principles to Clinical Applications in Drug Development

Abstract

This article provides a comprehensive overview of gene-based burden testing, a cornerstone method for identifying associations between rare genetic variants and complex traits. Tailored for researchers and drug development professionals, it covers foundational principles, core methodological approaches, and advanced optimization strategies for robust analysis. It further explores validation frameworks and comparative performance of leading tools like SAIGE, REGENIE, and SKAT-O, highlighting their application in large-scale biobanks for target discovery and safety assessment. The content synthesizes recent advances, including novel tests like SBAT and efficient meta-analysis methods, to equip scientists with the knowledge to implement these techniques effectively in genetic association studies.

Unraveling the Genetic Architecture: Why Rare Variants and Burden Tests Matter

Despite the success of genome-wide association studies (GWAS) in identifying thousands of common variants associated with complex diseases, a substantial portion of genetic heritability remains unexplained. Rare variants, typically defined as those with a minor allele frequency (MAF) below 0.5% to 1%, are increasingly recognized as crucial contributors to this "missing heritability" [1]. These variants often have larger effect sizes than common variants and are subject to purifying selection, making them particularly relevant for understanding disease pathogenesis [2]. Gene-based burden tests have emerged as powerful statistical frameworks for analyzing the collective contribution of multiple rare variants within a gene, overcoming power limitations inherent in single-variant analyses [1]. This protocol outlines comprehensive methodologies for identifying and validating the role of rare variants in complex diseases through advanced burden testing approaches.

Heritability Estimates and Study Findings

Table 1: Key Findings from Major Rare Variant Studies

Study / Finding	Sample Size	Disease Focus	Key Results
100,000 Genomes Project [3]	34,851 cases & family members	226 rare diseases	141 new disease-gene associations identified
Heritability Gap in T2D [1]	>150,000 individuals	Type 2 Diabetes	>70 GWAS loci explain only ~11% of heritability
Heritability Gap in Crohn's [1]	>210,000 individuals	Crohn's Disease	~70 GWAS loci explain only 23% of heritability
Hemiplegic Migraine Study [4]	184 patients	Hemiplegic Migraine	11 genes significantly associated via burden testing

Statistical Power Considerations for Rare Variant Analyses

Table 2: Statistical Methods for Rare Variant Analysis

Method Type	Examples	Best Use Case	Limitations
Burden Tests [2]	CAST, CMC, WSS	All variants in unit have similar effects	Power loss when both risk/protective variants exist
Variance-Component Tests [2]	SKAT, C-alpha	Mixed effect directions	Lower power when all variants have same effect
Combined Tests [2]	SKAT-O, Fisher	Unknown effect direction	Computational intensity
Case-Only Tests [5]	COBT	Limited control data	Relies on population mutation rates

Experimental Protocols for Gene-Based Burden Testing

Protocol 1: Case-Control Burden Analysis with geneBurdenRD

Application: Genome-wide rare variant association discovery in large cohorts [3]

Workflow:

Step-by-Step Methodology:

• Input Preparation

  • Collect VCF files from whole-genome or whole-exome sequencing of cases and controls
  • Assample metadata file with case-control assignments and covariates (age, sex, principal components)
  • Define analysis groups based on disease categories or phenotypic clusters [3]

• Variant Quality Control and Filtering

  • Remove samples with contamination indicators (unusually high heterozygosity)
  • Calculate quality metrics: read depth, transition/transversion ratio, novel variant counts
  • Apply variant-level filters: QUAL score, mapping quality, strand bias [1]
  • Retain variants with MAF < 0.5% in internal and external databases (gnomAD)

• Functional Annotation and Impact Prediction

  • Annotate variants using tools like Exomiser [3]
  • Predict functional impact with SIFT, PolyPhen, CADD for coding variants
  • Prioritize putative loss-of-function variants (nonsense, splice-site, frameshift)
  • Apply pathogenicity prediction algorithms for missense variants

• Variant Aggregation and Burden Calculation

  • Group variants by gene boundaries or functional units
  • Calculate burden metrics: variant counts, weighted sums based on predicted functionality
  • Adjust for covariates using regression frameworks

• Statistical Testing and Significance Assessment

  • Apply burden, variance-component, or omnibus tests as appropriate for hypothesis
  • Correct for multiple testing using Bonferroni or false discovery rate (FDR) methods
  • Validate associations in independent replication cohorts where possible

Protocol 2: Case-Only Burden Testing with COBT

Application: Rare disease research with limited control data [5]

Workflow:

Step-by-Step Methodology:

• Data Input and Processing

  • Input case-only sequencing data (whole exome or genome)
  • Obtain aggregated genotype data from public reference cohorts (gnomAD, 1000 Genomes)
  • Define rare variants based on population frequency thresholds (typically MAF < 0.1%)

• Variant Annotation and Filtering

  • Annotate variants with functional prediction scores
  • Filter based on quality metrics and sequencing depth
  • Apply inheritance model considerations (dominant, recessive, additive)

• Statistical Modeling with Poisson Framework

  • Model observed variant counts in cases using Poisson distribution
  • Calculate expected variant rates based on population mutation rates
  • Test for excess burden in cases versus population expectations
  • Account for multiple variants per individual and additive effects [5]

• Result Interpretation and Validation

  • Identify genes with significant excess burden of rare variants
  • Compare performance against alternative case-only methods
  • Validate known disease genes and prioritize novel candidates

Protocol 3: Multi-Omics Integration with Carrier Statistic

Application: Prioritizing functionally consequential rare variants by integrating expression data [6]

Workflow:

• Data Collection and Processing

  • Collect whole exome/genome sequencing data from cases and controls
  • Obtain matched transcriptomic data (RNA-seq) from relevant tissues
  • Perform standard QC on both genomic and transcriptomic datasets

• Variant Impact Quantification

  • Identify rare variants (MAF < 0.5%) in coding and regulatory regions
  • Quantify impact of rare variants on gene expression using association testing
  • Calculate carrier statistics measuring functional consequence magnitude

• Statistical Prioritization

  • Apply carrier statistic framework to prioritize variants with large functional effects
  • Test for enrichment of rare variants in cases versus controls
  • Integrate pathway information for biological interpretation

Application Notes and Case Studies

Case Study: 100,000 Genomes Project Application

The geneBurdenRD framework applied to 34,851 cases from the 100,000 Genomes Project identified 141 new disease-gene associations, with 69 prioritized after expert review [3]. Notable discoveries included:

• UNC13A in monogenic diabetes, implicating a known β-cell regulator
• GPR17 in schizophrenia, suggesting novel neuronal signaling mechanisms
• RBFOX3 in epilepsy, highlighting RNA splicing defects in neuronal excitability
• ARPC3 in Charcot-Marie-Tooth disease, connecting cytoskeletal regulation to neuropathy

This study demonstrated the clinical impact of large-scale statistical approaches, with potential to provide diagnoses for previously undiagnosed rare disease patients.

Case Study: Ciliopathy Research with COBT

Application of COBT to 478 ciliopathy patients successfully re-identified known causal genes and highlighted novel candidate genes in unsolved cases [5]. The case-only approach proved particularly valuable for rare diseases where matched controls are unavailable, effectively capturing additive effects and hypomorphic variants that alternative methods miss.

The Scientist's Toolkit: Essential Research Reagents

Table 3: Key Computational Tools and Resources for Rare Variant Burden Testing

Tool/Resource	Type	Primary Function	Application Context
geneBurdenRD [3]	R Package	Gene-based burden testing	Case-control studies in large cohorts
COBT [5]	Software Package	Case-only burden test	Studies lacking matched controls
Carrier Statistic [6]	Statistical Framework	Rare variant prioritization	Studies with multi-omics data
Exomiser [3]	Variant Prioritization	Functional annotation	Pre-processing for burden testing
gnomAD [5]	Reference Database	Population frequency data	Variant filtering and annotation
PathVar [4]	Bioinformatics Pipeline	Pathogenic variant identification	Variant prioritization in disease cohorts
1000 Genomes [5]	Reference Data	Population genetic variation	Control data for case-only analyses
Streptonigrin	Streptonigrin, CAS:1079893-79-0, MF:C25H22N4O8, MW:506.5 g/mol	Chemical Reagent	Bench Chemicals
TC-E 5008	6-Benzyl-1-Hydroxy-4-Methylpyridin-2(1H)-One	6-Benzyl-1-Hydroxy-4-Methylpyridin-2(1H)-One (SYC-435) is a potent, cell-active mutant IDH1 inhibitor for cancer research. This product is for Research Use Only (RUO) and is not intended for personal use.	Bench Chemicals

The investigation into the genetic underpinnings of human disease has undergone a fundamental transformation with the advent of next-generation sequencing (NGS). While genome-wide association studies (GWAS) successfully identified numerous common variants associated with complex diseases, they explained only a fraction of disease heritability and possessed limited power to detect rare variants [7]. This limitation catalyzed the development of gene-based collapsing methods, which aggregate rare variants from a genomic region of interest to test for association with diseases and traits [8]. The core premise of these methods is that multiple rare variants within a gene may collectively influence a shared trait, thereby amplifying signals that would be undetectable through single-variant analyses [9] [7].

These collapsing approaches have proven particularly valuable for identifying genes associated with both rare Mendelian disorders and complex diseases. For instance, in amyotrophic lateral sclerosis (ALS), collapsing analyses have pinpointed risk regions in known genes like SOD1, TARDBP, and FUS, while similar approaches have revealed novel cancer susceptibility genes including BIK, ATG12, and TG [9] [10]. The evolution from single-variant to gene-based analyses represents a critical methodological advancement, enabling researchers to systematically investigate the "rare variant-common disease" hypothesis across thousands of phenotypes in large biobanks [11].

Fundamental Principles of Gene-Based Collapsing Analyses

Core Concept and Rationale

Gene-based collapsing analyses operate on a fundamental principle: instead of testing individual rare variants for association with a disease, these methods aggregate multiple rare variants within a predefined genomic unit—typically a gene or functional domain—and test whether the burden of these variants is significantly higher in cases than in controls [9] [8]. This approach addresses the critical power limitation of single-variant tests for rare variants, which may be observed too infrequently to demonstrate statistical significance individually [8] [7].

The biological rationale underpinning these methods stems from the observation that pathogenic mutations in many diseases are concentrated in specific genic regions, particularly functional domains critical to protein function [9]. Furthermore, rare variants with potentially large effect sizes are more likely to be subject to purifying selection, making them particularly relevant for disease investigation [8]. By aggregating these variants, collapsing methods can detect genes where multiple different rare variants contribute to disease risk, even when no single variant reaches significance.

Key Methodological Variations

Several methodological frameworks have been developed for gene-based collapsing analyses, each with distinct approaches to variant aggregation and statistical testing:

Table 1: Categories of Gene-Based Collapsing Methods

Method Category	Key Principle	Advantages	Limitations
Burden Tests	Aggregates variants into a single score assuming unidirectional effects [8] [7]	High power when most variants are causal with same effect direction	Power loss with non-causal variants or bidirectional effects [7]
Variance Component Tests (e.g., SKAT)	Tests distribution of variant effects without assuming unidirectional effects [8] [7]	Robust to bidirectional effects and presence of non-causal variants	Lower power when most variants are causal with same direction [7]
Combined Tests (e.g., SKAT-O)	Hybrid approach combining burden and variance component tests [7]	Balanced performance across different genetic architectures	Computational complexity [7]
Domain-Based Collapsing	Focuses burden tests on functional protein domains rather than entire genes [9]	Increased power for genes with regionally localized pathogenic variation	Requires prior knowledge of functional domains

Experimental Design and Workflow

Defining Qualifying Variants

A critical step in any collapsing analysis is establishing criteria for "qualifying variants" to include in the aggregate test. This process typically involves multiple filtering layers based on functional annotation, population frequency, and predicted functional impact:

• Functional Annotation: Variants are typically restricted to those likely to impact protein function, including loss-of-function (LoF) variants (nonsense, frameshift, canonical splice site), and often missense variants with predicted deleterious effects using tools like PolyPhen-2 or SIFT [9] [11].

• Allele Frequency Threshold: Variants are filtered based on maximum allele frequency, typically with minor allele frequency (MAF) < 0.1% to 1% in population databases like gnomAD, to focus on rare variation [9] [11] [12].

• Quality Control: Rigorous variant quality control is essential, including filters for call rate, sequencing depth, and genotype quality to minimize false positives [9] [3].

Different studies employ various combinations of these criteria. For example, a study of amyotrophic lateral sclerosis defined qualifying variants as "nonsynonymous coding or canonical splice variants that have a minor allele frequency (MAF) ≤0.1% in cases and controls" [9], while cancer susceptibility research has aggregated "rare missense and loss of function variants" [10].

Analytical Workflow

The following diagram illustrates the standard workflow for gene-based collapsing analysis:

Statistical Considerations and Multiple Testing

Appropriate statistical handling is crucial for valid collapsing analyses. Several key considerations include:

• Handling Population Stratification: Methods like principal component analysis (PCA) are essential to correct for population stratification that can cause spurious associations [7]. Linear mixed models (LMMs) can also account for relatedness and population structure [11].

• Case-Control Imbalance: For binary traits with severe case-control imbalance, applying a minimum expected carrier threshold (e.g., at least 10 variant carriers in cases) helps avoid false positives [11].

• Multiple Testing Correction: Given the testing of thousands of genes, stringent significance thresholds are required. Study-wide significance levels of ( p < 3.4 × 10^{-10} ) or gene-based Bonferroni correction (( p < 2.5 × 10^{-6} ) for 20,000 genes) are commonly applied [11] [12].

Advanced Applications and Methodological Innovations

Domain-Based Collapsing Approaches

Standard gene-based collapsing methods that treat all variants across a gene as equivalent can miss signals when pathogenic mutations cluster in specific functional regions. Domain-based collapsing addresses this limitation by using homology-defined protein domains as the unit for collapsing analysis rather than entire genes [9].

This approach has proven particularly powerful for genes where pathogenic variants are concentrated in specific domains. In ALS research, while standard gene collapsing identified SOD1, it failed to detect significant associations for TARDBP and FUS. However, domain-based collapsing revealed genome-wide significant associations for a glycine-rich domain in TARDBP (OR = 7; P = 5.84 × 10⁻⁷) and an Arg-Gly rich domain in FUS (OR = 8.6; P = 3.6 × 10⁻⁵), highlighting the increased resolution of this approach for genes with regionally localized pathogenic variation [9].

Case-Only Study Designs

The challenge of obtaining matched controls in rare disease research has spurred the development of case-only burden tests. The recently developed Case-Only Burden Test (COBT) uses a Poisson model to test for excess variants in a gene compared to expectations from general population mutation rates, effectively using public reference cohorts as controls [5].

Applied to 478 ciliopathy patients, COBT successfully re-identified known causal genes and highlighted novel candidate variants in unsolved cases, demonstrating that case-only designs can enable gene discovery when traditional case-control analyses are not feasible [5]. This approach is particularly valuable for rare diseases where recruiting matched controls is challenging.

Large-Scale Applications in Biobanks

The integration of collapsing analyses with large-scale biobanks has dramatically expanded the scope of rare variant research. One comprehensive analysis applied gene-based collapsing to 4,264 phenotypes in 49,960 exome-sequenced individuals from the UK Biobank, complemented by 1,934 phenotypes in 21,866 participants from the Healthy Nevada Project [11].

This study identified 64 statistically significant gene-based associations in the meta-analysis, demonstrating that singletons (variants observed only once) made significant contributions to the results. Importantly, the vast majority of these associations could not have been identified with genotyping chip data alone, highlighting the unique power of sequencing-based rare variant detection [11].

Table 2: Key Research Reagents and Computational Tools for Gene-Based Collapsing Analyses

Tool/Resource	Type	Primary Function	Application Notes
Exomiser [3]	Software	Variant prioritization tool	Used in the 100,000 Genomes Project to identify rare putative disease-causing variants for burden testing
geneBurdenRD [3]	R Framework	Gene burden testing for rare diseases	Open-source R framework specifically designed for rare disease sequencing cohorts
COBT [5]	Software	Case-Only Burden Test	Implements Poisson model for case-only designs using public reference cohorts as controls; available at https://github.com/RausellLab/COBT
collapseRows [13]	R Function	Data aggregation and collapsing	Implements standard and network-based methods for collapsing multiple related variables into representatives
SKAT/SKAT-O [8] [7]	Statistical Package	Rare variant association testing	Robust methods that accommodate bidirectional variant effects and presence of non-causal variants
gnomAD [5] [11]	Database	Population variant frequencies	Essential for filtering out common variants and establishing frequency thresholds for qualifying variants
Conserved Domain Database (CDD) [9]	Database	Protein domain definitions	Used for domain-based collapsing approaches to define functional regions within genes

Integrated Workflow for Comprehensive Gene Burden Analysis

The following diagram illustrates an integrated workflow for comprehensive gene burden analysis, incorporating both standard and advanced approaches:

Protocol: Implementing a Gene-Based Collapsing Analysis

Sample Preparation and Sequencing

• Cohort Selection: Define cases and controls based on precise phenotypic criteria. For rare diseases, consider including family members to enhance power [3].
• DNA Extraction: Use high-quality DNA extraction protocols suitable for whole exome or genome sequencing.
• Library Preparation and Sequencing: Perform whole exome capture or whole genome sequencing using established platforms (Illumina, for example) to achieve sufficient coverage (typically >30x for WES, >15x for WGS).

Variant Calling and Quality Control

• Variant Calling: Process raw sequencing data through standard pipelines (BWA-GATK, for example) to generate variant call format (VCF) files.
• Quality Filtering: Apply stringent quality filters including call rate >95%, genotype quality >20, and read depth >10 to minimize false positives.
• Sample QC: Remove samples with excessive heterozygosity, contamination, or mismatched phenotypic-genetic sex.
• Variant Annotation: Annotate variants using tools like ANNOVAR or VEP with functional predictions (SIFT, PolyPhen-2), population frequency databases (gnomAD), and transcript information.

Gene-Based Collapsing Analysis

• Define Qualifying Variants: Establish criteria based on:

  • Minor Allele Frequency: Typically MAF < 0.1% in population databases [11]
  • Functional Impact: Loss-of-function (stop-gain, frameshift, canonical splice site) and deleterious missense variants [10] [11]
  • Quality Metrics: Ensure high-quality variant calls

• Perform Association Testing:

  • For burden tests: Aggregate qualified variants per gene and test for frequency differences between cases and controls using regression models adjusted for covariates [8] [7]
  • For SKAT/SKAT-O: Test for associations using variance-component models that accommodate bidirectional effects [8] [7]
  • For domain-based analyses: Restrict collapsing to specific functional domains defined by databases like CDD [9]

• Address Population Stratification:

  • Calculate principal components (PCs) from common variants
  • Include top PCs as covariates in association models [7]
  • Consider using linear mixed models to account for relatedness [11]

Significance Assessment and Validation

• Multiple Testing Correction: Apply Bonferroni correction based on the number of genes tested or false discovery rate (FDR) control [12]
• Replication Analysis: Seek independent replication in additional cohorts where possible [11]
• Functional Validation: Plan experimental follow-up using model systems or functional assays for novel associations

Interpretation and Clinical Translation

Successful application of gene-based collapsing analyses has yielded significant insights across diverse diseases. In the 100,000 Genomes Project, this approach identified 141 new disease-gene associations for rare Mendelian diseases, with 30 associations supported by existing experimental evidence [3]. Similarly, cancer burden analyses have revealed both risk-increasing genes (BIK for prostate cancer) and protective genes (PPP1R15A for breast cancer), suggesting potential therapeutic targets [10].

When interpreting results, researchers should consider that different collapsing methods may yield complementary insights. For example, in ALS research, a standard gene-based approach identified SOD1, while domain-based methods were necessary to detect significant associations for TARDBP and FUS [9]. Similarly, case-only designs like COBT offer viable alternatives when matched controls are unavailable [5].

As sequencing costs continue to decline and sample sizes grow, gene-based collapsing analyses will play an increasingly central role in unraveling the genetic architecture of both rare and common diseases, ultimately accelerating therapeutic development and improving diagnostic yield for patients with undiagnosed genetic conditions.

Gene-based burden tests represent a cornerstone of modern rare variant analysis, addressing a critical limitation of single-variant approaches by aggregating the effects of multiple rare genetic variants within a defined gene or genomic region. The fundamental premise underlying these tests is that complex phenotypes may arise not from single deleterious variants, but from the collective burden of multiple rare variants within biologically relevant units. This approach has proven particularly valuable in scenarios where individual rare variants occur too infrequently to achieve statistical significance in association testing, yet cumulatively contribute to disease pathogenesis through haploinsufficiency, dominant-negative effects, or other mechanisms. The evolution of burden testing methodologies has paralleled advances in sequencing technologies and the availability of large-scale genomic resources, enabling applications across the spectrum of human disease from monogenic disorders to complex polygenic traits.

The statistical foundation of burden tests involves comparing the aggregate frequency of rare variants in cases versus controls, operating under the assumption that affected individuals carry more deleterious variants in disease-relevant genes. Early implementations simply counted the number of rare variants per gene, while contemporary approaches incorporate functional annotations, allele frequency thresholds, and directional effects to improve power and specificity. As genomic medicine progresses, burden tests have expanded beyond research settings into diagnostic pipelines, where they facilitate the interpretation of rare variants in patients with undiagnosed diseases and help prioritize candidate genes for functional validation.

Application in Mendelian Disorders

Overcoming Challenges in Rare Disease Genetics

Mendelian disorders present unique challenges for genetic analysis, including extensive locus heterogeneity, incomplete penetrance, and the scarcity of affected individuals within individual studies. Gene-based burden tests have emerged as powerful tools to address these limitations by enabling the aggregation of signals across multiple variants and individuals. Traditional linkage analysis and single-variant association tests often fail in genetically heterogeneous conditions or when dealing with de novo mutations, making burden approaches particularly valuable for these scenarios.

Table 1: Key Studies Applying Burden Tests to Mendelian Disorders

Study	Disease Focus	Sample Size	Key Findings	Methodological Innovations
COBT [5]	Ciliopathies	478 patients	Re-identified known causal genes and highlighted novel candidate variants	Case-only design using public reference cohorts; Poisson model accounting for multiple variants per individual
TRAPD [14]	Idiopathic Hypogonadotropic Hypogonadism (IHH)	393 cases vs. 123,136 gnomAD controls	Rediscovered known IHH genes (FGFR1, TACR3, GNRHR); identified novel association with TYRO3	Adaptive variant quality filtering; use of synonymous variants for calibration
Exomiser Optimization [15]	Various rare diseases (UDN cohort)	386 diagnosed probands	Improved diagnostic variant ranking in top 10 from 49.7% to 85.5% for GS data	Parameter optimization for variant prioritization tools

In idiopathic hypogonadotropic hypogonadism (IHH), a rare Mendelian disorder with significant locus heterogeneity, researchers demonstrated the utility of burden testing by analyzing whole-exome sequencing data from 393 cases against 123,136 controls from gnomAD [14]. This approach successfully re-discovered known IHH genes including FGFR1, TACR3, and GNRHR, while also identifying a significant burden in TYRO3, a gene previously implicated in hypogonadotropic hypogonadism in mouse models. The study introduced the TRAPD software, which implements sophisticated variant quality filtering to address technical artifacts when using public control databases.

Case-Only Study Designs

For many rare diseases, recruiting matched controls presents practical and financial challenges. The recently developed Case-Only Burden Test (COBT) addresses this limitation by using aggregated genotypes from public reference cohorts, employing a Poisson model to test for excess variants in a gene compared to expectations based on general population mutation rates [5]. When applied to 478 ciliopathy patients, COBT successfully re-identified known causal genes and highlighted novel candidate variants in previously unsolved cases, demonstrating sensitivity comparable to case-control designs while overcoming control recruitment barriers.

The diagnostic utility of burden testing is further enhanced when integrated with multi-omic approaches. Recent research from Baylor Genetics highlights how RNA sequencing provides functional evidence to reclassify variants of uncertain significance identified through exome and genome sequencing [16]. In their cohort of 3,594 consecutive cases, RNA-seq enabled reclassification of half of eligible variants, with over one-third of RNA-seq eligible cases having noncoding variants that would have been missed by exome sequencing alone.

Variant Prioritization in Diagnostic Pipelines

In clinical diagnostics, burden tests have been integrated into variant prioritization pipelines to improve diagnostic yield. A systematic evaluation of the Exomiser/Genomiser software suite on 386 diagnosed probands from the Undiagnosed Diseases Network demonstrated that parameter optimization could dramatically improve performance, increasing the percentage of coding diagnostic variants ranked within the top 10 candidates from 49.7% to 85.5% for genome sequencing data, and from 67.3% to 88.2% for exome sequencing data [15]. This optimized approach has been implemented in the Mosaic platform to support ongoing analysis of undiagnosed participants, providing an efficient, scalable framework for periodic reanalysis.

Application in Complex Disease Genomics

Expanding to Polygenic and Common Diseases

While initially developed for Mendelian disorders, gene-based burden tests have increasingly been applied to complex diseases, where they capture the contribution of rare variants with potentially larger effect sizes than common variants identified through genome-wide association studies (GWAS). The underlying hypothesis is that rare coding variants with moderate to large effects contribute to the "missing heritability" not explained by common variants.

Table 2: Burden Test Applications in Complex Diseases and Biomarkers

Study	Phenotype Category	Sample Size	Key Genes Identified	Effect Directions
Cancer Susceptibility [10]	22 cancer sites	130,991 cases vs. 733,486 controls	BIK (prostate), ATG12 (colorectal), TG (thyroid), CMTR2 (lung/melanoma)	Increased risk
			AURKB (any cancer), PPP1R15A (breast)	Decreased risk
Blood Biomarkers [17]	28 blood biomarkers	UK Biobank cohort	ALPL (alkaline phosphatase), LDLR (LDL), PCSK9 (LDL)	Positive and negative effects observed
Meta-SAIGE [18]	83 low-prevalence phenotypes	UK Biobank + All of Us	237 gene-trait associations	Mixed

In cancer genomics, a large-scale burden analysis across 22 cancer sites involving 130,991 cases and 733,486 controls from Iceland, Norway, and the United Kingdom identified six cancer susceptibility genes [10]. The study found four genes associated with increased cancer risk: BIK (prostate cancer), ATG12 (colorectal cancer), TG (thyroid cancer), and CMTR2 (both lung cancer and cutaneous melanoma). Additionally, two genes showed protective effects, with rare variants in AURKB associated with decreased risk for any cancer and PPP1R15A variants associated with reduced breast cancer risk. These findings highlight biological pathways involving apoptosis, autophagy, and cell stress response as potential targets for therapeutic development.

Integration with Common Variant Approaches

Complex traits are influenced by both common and rare variants, yet these are typically analyzed independently. Recent research has explored integrated models that combine polygenic risk scores (PRS) based on common variants with gene-based burden scores from rare variants [17]. When applied to 28 blood biomarkers in the UK Biobank dataset, association analyses revealed genes with significant effects, including ALPL for alkaline phosphatase (z-score: -49.6) and LDLR for LDL direct measurement (z-score: 23.4). The directionality of these effects provided insights into biological mechanisms: damaging variants in ALPL decreased alkaline phosphatase levels, while damaging variants in LDLR increased LDL levels, consistent with familial hypercholesterolemia pathogenesis.

Interestingly, despite strong individual-level effects of rare variants, combined prediction models incorporating both PRS and gene-based scores showed limited improvement over PRS-only models for most biomarkers, suggesting that common variant-based PRS might be more informative for population-level risk prediction [17]. This highlights the complementary nature of common and rare variants: while rare variants can have dramatic effects at the individual level, their population-level contribution to complex trait variability may be limited by their low frequency.

Experimental Protocols and Methodologies

Standard Burden Test Protocol

The following protocol outlines a standardized approach for conducting gene-based burden tests, synthesizing methodologies from multiple studies [5] [14] [17]:

Step 1: Variant Calling and Quality Control

• Perform whole-exome or whole-genome sequencing using standard platforms (Illumina recommended)
• Align sequences to reference genome (GRCh38 preferred) using BWA-MEM
• Conduct variant calling following GATK Best Practices, including base quality score recalibration, indel realignment, and duplicate removal
• Implement strict quality filters: read depth >10×, genotype quality >20, and missingness <5%
• Annotate variants using Ensembl VEP with additional pathogenicity predictors (PolyPhen-2, SIFT, CADD)

Step 2: Variant Selection and Filtering

• Define rare variant threshold (typically MAF < 0.1% in gnomAD or population-matched databases)
• Filter for protein-altering variants: loss-of-function (stop-gain, frameshift, splice-site), missense, and in-frame indels
• Exclude variants in low-complexity regions and segmental duplications
• Apply functional impact predictions to prioritize deleterious variants

Step 3: Gene-Based Aggregation

• Group qualifying variants by gene boundaries (GENCODE definitions)
• Calculate burden scores for each individual: simple count, weighted by functional impact, or allele count
• Adjust for gene-specific mutation rates and coverage

Step 4: Association Testing

• For case-control designs: Use logistic regression with burden score as predictor, adjusting for covariates (principal components, sequencing platform)
• For case-only designs: Compare observed variant counts to expected counts based on population mutation rates [5]
• Implement significance thresholds adjusted for multiple testing (Bonferroni correction for gene count or false discovery rate)

Step 5: Validation and Replication

• Replicate findings in independent cohorts when available
• Conduct sensitivity analyses with different variant inclusion criteria
• Perform functional validation through experimental approaches or orthogonal datasets

Meta-Analysis Protocol for Multi-Cohort Studies

Meta-analysis substantially enhances power for rare variant association tests by combining summary statistics across cohorts. The Meta-SAIGE protocol provides a scalable framework [18]:

Step 1: Cohort-Level Preparation

• Generate per-variant score statistics and sparse LD matrices for each cohort using SAIGE
• Adjust for case-control imbalance using saddlepoint approximation
• Account for sample relatedness using genetic relationship matrices

Step 2: Summary Statistics Harmonization

• Combine score statistics from all studies into a single superset
• Recalculate variances by inverting SPA-adjusted p-values
• Apply genotype-count-based SPA for improved type I error control

Step 3: Gene-Based Testing

• Conduct Burden, SKAT, and SKAT-O tests on harmonized summary statistics
• Collapse ultrarare variants (MAC < 10) to enhance power
• Combine p-values across functional annotations and MAF cutoffs using Cauchy combination method

This approach has been shown to effectively control type I error rates even for low-prevalence binary traits (1-5%) and achieves power comparable to pooled individual-level data analysis [18].

Signaling Pathways and Biological Mechanisms

Gene-based burden tests have illuminated key biological pathways contributing to human diseases through the aggregation of rare variant signals:

Figure 1: Biological Pathways Implicated by Burden Tests. Rare variant burden in key genes affects disease risk through disruption of fundamental cellular processes.

The pathways highlighted through burden analyses reveal convergent mechanisms across seemingly disparate conditions. In cancer, autophagy (ATG12), apoptosis (BIK), and stress response (PPP1R15A) pathways emerge as critical protective mechanisms [10]. For endocrine disorders like IHH, hormonal signaling pathways (GNRHR, TACR3) represent key vulnerability points [14]. In metabolic traits, cholesterol homeostasis genes (LDLR, PCSK9) demonstrate how rare variants with large effects can inform both pathophysiology and therapeutic targeting [17].

The Scientist's Toolkit: Essential Research Reagents

Table 3: Key Research Reagents and Computational Tools for Burden Testing

Tool/Resource	Type	Primary Function	Application Context
COBT [5]	Software	Case-only burden test using public references	Mendelian disorders with limited controls
TRAPD [14]	Software	Burden testing with public database controls	Rare disease gene discovery
Meta-SAIGE [18]	Software	Rare variant meta-analysis	Multi-cohort studies of complex traits
Exomiser/Genomiser [15]	Software	Variant prioritization integrating phenotype	Diagnostic pipeline optimization
gnomAD	Database	Population variant frequencies	Variant filtering and annotation
UK Biobank	Resource	Integrated genotype and phenotype data	Complex trait analyses
VEP	Tool	Variant effect prediction	Functional annotation of variants
SAIGE-GENE+ [18]	Software	Rare variant tests for biobank data	Large-scale association testing
Spantide I	Spantide I, MF:C75H108N20O13, MW:1497.8 g/mol	Chemical Reagent	Bench Chemicals
Valone	Valone, CAS:83-28-3, MF:C14H14O3, MW:230.26 g/mol	Chemical Reagent	Bench Chemicals

Successful implementation of gene-based burden tests requires careful selection of computational tools and biological resources. The field has evolved from simple counting methods to sophisticated frameworks that address specific study design challenges. For case-only designs, COBT provides a specialized solution that leverages public reference data [5], while TRAPD facilitates the use of gnomAD as control data with appropriate quality controls [14]. For large-scale collaborations, Meta-SAIGE enables meta-analysis while controlling type I error and maintaining computational efficiency [18].

In diagnostic settings, Exomiser and Genomiser represent critical tools for variant prioritization, with optimized parameters significantly improving performance [15]. The integration of phenotypic information through Human Phenotype Ontology terms further enhances the specificity of these approaches for rare disease diagnosis. For complex trait analysis, methods that combine burden scores with polygenic risk scores enable a more comprehensive assessment of genetic architecture [17].

Workflow and Experimental Design

The implementation of gene-based burden tests follows a structured workflow that ensures robust and interpretable results:

Figure 2: Gene-Based Burden Test Workflow. Key considerations at each step ensure robust results and minimize technical artifacts.

Critical decision points throughout this workflow significantly impact the results and interpretation of burden tests. Population stratification represents a particular challenge, requiring careful adjustment using principal components or genetic relationship matrices [14]. Coverage uniformity across genes and samples must be addressed, potentially through depth-based weighting or exclusion of poorly covered regions. The selection of functional annotation criteria and MAF thresholds should be guided by the specific disease model, with more restrictive filters increasing specificity but potentially reducing sensitivity for genes tolerant to certain mutation types.

Multiple testing correction remains challenging due to the correlation between tests (genes are not independent) and the implementation of multiple burden tests with different parameters. Approaches such as Bonferroni correction, false discovery rate control, or permutation testing provide different tradeoffs between stringency and power. Finally, independent replication and functional validation remain essential for confirming associations, particularly for novel gene-disease relationships.

The fundamental statistical challenge in rare variant association analysis stems from the low minor allele frequency (MAF) of the variants under investigation. For single genetic variants that are rare (conventionally defined as MAF < 0.5% to 1%), the statistical power to detect their individual association with a trait using traditional genome-wide association study (GWAS) methods is severely limited, even in very large sample sizes [19] [20]. This limitation arises because the expected number of observed minor alleles in a sample is low, making it difficult to distinguish true associations from background variation with statistical confidence. Consequently, while common variants typically require single-variant tests with massive sample sizes to achieve genome-wide significance, this approach becomes practically infeasible for rare variants due to prohibitive sample size requirements that often exceed millions of individuals [19].

This statistical reality has motivated the development of gene-based aggregation methods that collectively test the association of multiple rare variants within biologically meaningful units, most commonly genes [21] [20]. The core rationale behind these methods is that while individual rare variants may be too rare to detect individually, their collective impact within functional genomic regions (such as genes or pathways) may yield a detectable signal when analyzed jointly. This approach aligns with the biological hypothesis that multiple rare variants within the same gene may disrupt its function through different molecular mechanisms but ultimately contribute to similar phenotypic outcomes [19] [7]. Gene-based burden testing has proven particularly successful in identifying novel cancer susceptibility genes and rare disease associations in large-scale sequencing studies [10] [3].

Table 1: Key Challenges in Rare Variant Association Studies

Challenge	Impact on Statistical Power	Potential Solutions
Low allele frequency	Reduced power for single-variant tests	Gene-based aggregation methods
Allelic heterogeneity	Multiple causal variants with different effect sizes	Burden tests, SKAT, SKAT-O
Presence of non-causal variants	Noise in association signal	Functional annotation filtering
Bidirectional effects	Signal cancellation in burden tests	Variance-component tests (SKAT)
Case-control imbalance	Type I error inflation	Saddlepoint approximation methods

Methodological Approaches for Overcoming Power Limitations

Burden Tests: The Collapsing Principle

Burden tests represent the foundational approach to rare variant aggregation, operating on the principle of collapsing multiple variants into a single genetic score for each individual [20]. The core methodology involves creating a burden score that summarizes the cumulative rare variant load within a predefined genomic region (typically a gene) for each study participant. This score can be a simple count of the number of rare alleles carried (assuming an additive model) or a weighted sum that incorporates variant-specific weights, often based on functional predictions or frequency [7] [20]. The statistical association between this burden score and the phenotype of interest is then tested using regression frameworks, with the null hypothesis being that the collective burden of rare variants in the region is not associated with the trait [20].

The most straightforward burden test is the Cohort Allelic Sums Test (CAST), which collapses variants into a binary indicator of carrier status (presence or absence of any qualifying rare variant) and tests for association using Fisher's exact test or logistic regression [22] [20]. More sophisticated approaches include the Weighted Sum Statistic (WSS), which incorporates weights inversely related to variant frequency, giving more weight to rarer variants under the assumption that they may have larger effect sizes [22]. While burden tests are powerful when most aggregated variants are causal and influence the trait in the same direction, their power diminishes substantially when these assumptions are violated, particularly in the presence of non-causal variants or variants with bidirectional effects (both risk and protective) [21] [23].

Variance-Component Tests: Accommodating Effect Heterogeneity

Variance-component tests, most notably the Sequence Kernel Association Test (SKAT), address a key limitation of burden tests by allowing for heterogeneous effect directions and sizes among the aggregated rare variants [18] [20]. Instead of collapsing variants into a single score, SKAT models the variant-phenotype relationship using a random effects framework where each variant is assumed to have its own effect size, following a distribution with a mean of zero and a common variance [22] [20]. The test statistic is based on a weighted sum of the squares of the scores for individual variants, making it robust to the presence of both non-causal variants and variants with opposing effects on the trait [22].

SKAT is particularly advantageous in scenarios where the genetic architecture of a trait includes both risk and protective variants within the same gene, or when only a small proportion of the aggregated variants are truly causal [21] [23]. The method employs a kernel function to model the genetic similarity between individuals based on their rare variant profiles, and tests whether this similarity correlates with phenotypic similarity [20]. Empirical studies have demonstrated that SKAT maintains better power than burden tests when the proportion of causal variants is low (<30%) or when variants have bidirectional effects, though it may be less powerful when most variants are causal and have effects in the same direction [21] [23].

Adaptive and Omnibus Tests: Strategic Hybridization

Recognizing that the optimal statistical approach depends on the underlying genetic architecture—which is typically unknown a priori—researchers have developed adaptive methods that combine the strengths of burden and variance-component approaches [18] [20]. The most widely used omnibus test is SKAT-O, which optimally combines burden and SKAT statistics using a data-driven approach that estimates the correlation between variant effects [22] [20]. SKAT-O includes a tuning parameter, ρ, that ranges from 0 to 1, where ρ = 0 corresponds to the standard SKAT test and ρ = 1 corresponds to a burden test [20]. The method estimates this parameter from the data to maximize power, effectively adapting to the true genetic architecture.

Another adaptive approach is the Aggregated Cauchy Association Test (ACAT), which combines p-values from different aggregation methods or from tests with different variant weightings [7]. This method is particularly useful when integrating evidence across multiple functional annotations or MAF thresholds, as it does not require assumptions about the correlation structure between tests and remains powerful when only a small number of variants are causal [7]. These adaptive methods have demonstrated robust performance across diverse genetic architectures, making them particularly valuable for hypothesis-free scanning where the true architecture is unknown [21] [20].

Table 2: Comparison of Major Rare Variant Association Methods

Method	Class	Key Assumptions	Optimal Use Case
CAST	Burden	All variants causal, same effect direction	High proportion of causal variants with uniform effects
Weighted Sum	Burden	Rare variants have larger effects	Causal variants with effect sizes inversely correlated to frequency
SKAT	Variance-component	Heterogeneous effect directions	Mixed protective/risk variants or low proportion of causal variants
SKAT-O	Omnibus	Adapts to underlying architecture	Unknown genetic architecture or mixed scenarios
ACAT	P-value combination	Combines complementary signals	Integration across multiple annotations or MAF thresholds

Practical Implementation and Protocol Guidance

Sample Size and Study Design Considerations

The statistical power of rare variant association studies is profoundly influenced by both total sample size and the balance between cases and controls, particularly for binary traits [23]. Empirical power simulations have revealed that for balanced case-control designs (1:1 ratio), SKAT generally achieves higher power than burden tests with sample sizes exceeding 4,000 individuals, whereas burden tests may perform better in smaller sample sizes (<1,000 total) when the proportion of causal variants is high [23]. For unbalanced designs, which are common for rare diseases, SKAT typically demonstrates superior performance, with simulations indicating that case numbers rather than case-control ratios primarily drive power when control groups are large [23].

Notably, achieving 90% power to detect rare variant associations with moderate effect sizes (OR ≈ 2.5) generally requires at least 200 cases for SKAT in unbalanced designs with large control groups (e.g., 10,000 controls), whereas burden tests may require 500-1,000 cases to achieve comparable power under similar conditions [23]. These sample size requirements are further influenced by the proportion of causal variants aggregated in the test, with burden tests requiring a higher proportion of causal variants (>50-60%) to maintain advantage over variance-component methods [21]. For very rare variants (MAF < 0.01%), even larger sample sizes are typically necessary, motivating the increasing emphasis on meta-analytic approaches across multiple biobanks [18].

Variant Filtering and Annotation Strategies

The power of gene-based rare variant tests can be substantially improved through strategic variant filtering and the incorporation of functional annotations [7] [3]. A critical first step involves defining the "rare" threshold, with current conventions typically using MAF cutoffs of 0.1-1% for complex traits and even lower thresholds (0.01-0.05%) for Mendelian diseases [20]. Beyond frequency filtering, focusing on putative functional categories—such as protein-truncating variants (PTVs), deleterious missense variants, or variants in specific functional domains—can enrich for causal variants and improve power [21] [3].

The implementation of variant masks that specify which categories of rare variants to include in aggregation tests has been shown to significantly impact power [21]. For example, restricting analyses to PTVs and deleterious missense variants (as defined by computational prediction algorithms) typically increases power when these variant classes have higher probabilities of being causal [21]. Additionally, incorporating functional annotations and intolerance metrics (e.g., missense constraint Z-scores) as weights in burden tests or SKAT can further improve power by upweighting variants more likely to have functional consequences [7] [3].

Figure 1: Rare Variant Analysis Workflow

Statistical Analysis Protocol for Gene-Based Burden Testing

A comprehensive protocol for gene-based burden analysis involves multiple stages of data processing and statistical testing. The following protocol outlines key steps for conducting robust rare variant association studies:

Step 1: Quality Control and Variant Calling Begin with high-quality whole exome or genome sequencing data processed through standardized pipelines. Implement quality control filters to remove samples with low call rates, contamination, or discordant phenotypic information. Apply variant-level filters to remove low-quality calls using metrics such as depth, genotype quality, and missingness [3] [24]. For population structure control, generate principal components from common variants and include them as covariates in association models [7].

Step 2: Variant Annotation and Functional Filtering Annotate remaining variants using functional prediction algorithms (e.g., SIFT, PolyPhen, CADD) and population frequency databases (e.g., gnomAD). Define variant masks based on predicted functional impact—typical masks include protein-truncating variants only, PTVs plus deleterious missense, or all rare missense variants [21] [3]. Consider gene-level constraint metrics (e.g., pLI scores) to prioritize genes intolerant to variation.

Step 3: Burden Score Calculation For each sample and genomic region (typically gene-based), calculate burden scores according to the chosen masking scheme. Common approaches include: (1) binary indicators of carrier status for any qualifying variant; (2) weighted counts of rare alleles with weights based on functional prediction scores or frequency; or (3) non-synonymous variant counts stratified by predicted functional impact [10] [3].

Step 4: Association Testing Test the association between burden scores and phenotypes using appropriate statistical models. For continuous traits, employ linear regression; for binary traits, use logistic regression with Firth's correction or saddlepoint approximation to address case-control imbalance [18] [23]. Include relevant covariates such as age, sex, sequencing platform, and genetic principal components to control for confounding.

Step 5: Significance Evaluation and Multiple Testing Correction Assess statistical significance using gene-based thresholds, typically employing Bonferroni correction based on the number of genes tested (e.g., exome-wide significance threshold of 2.5×10^-6 for 20,000 genes) [10] [3]. For hypothesis-driven analyses with predefined gene sets, less stringent thresholds may be appropriate.

Advanced Methodological Considerations

Addressing Population Stratification and Technical Artifacts

Population stratification poses particular challenges for rare variant association studies, as rare variants often have geographically restricted distributions and can create spurious associations if ancestry is correlated with both variant frequency and phenotype [19] [20]. Standard approaches for controlling stratification in common variant analyses, such as genomic control and principal component adjustment, may be insufficient for rare variants due to their low frequencies and limited contributions to ancestry inference [20]. More sophisticated methods, including mixed models that incorporate a genetic relatedness matrix (GRM) and methods that specifically model rare variant population structure, are often necessary to adequately control type I error [18] [20].

Technical artifacts from sequencing or variant calling can also disproportionately impact rare variant analyses, as errors may be confused with genuine rare variants [19]. Differential missingness or quality between cases and controls can introduce spurious associations, necessitating rigorous quality control and potentially the use of methods that account for uncertainty in variant calling [3]. For binary traits with extreme case-control imbalances, standard asymptotic tests (e.g., score tests) can exhibit substantial type I error inflation, requiring specialized approaches such as saddlepoint approximation (implemented in SAIGE and Meta-SAIGE) or Firth's correction to maintain proper error control [18] [23].

Meta-Analysis Methods for Collaborative Studies

As single cohorts rarely provide sufficient sample sizes for well-powered rare variant analyses, meta-analysis methods have become essential for combining evidence across multiple studies [18]. Recent advances in rare variant meta-analysis include methods such as Meta-SAIGE, which extends the SAIGE framework to summary-level data while maintaining accurate type I error control even for low-prevalence binary traits [18]. Unlike earlier approaches like MetaSTAAR, Meta-SAIGE implements a two-level saddlepoint approximation that effectively addresses case-control imbalance and allows reuse of linkage disequilibrium matrices across phenotypes, significantly improving computational efficiency for phenome-wide analyses [18].

Empirical evaluations have demonstrated that meta-analysis approaches can identify numerous associations that are not detectable in individual cohorts alone. For example, a Meta-SAIGE analysis of 83 low-prevalence phenotypes in UK Biobank and All of Us data identified 237 gene-trait associations, 80 of which were not significant in either dataset individually [18]. These findings highlight the critical importance of collaborative efforts and sophisticated meta-analysis methods for comprehensively characterizing the contribution of rare variants to human traits and diseases.

Figure 2: Rare Variant Meta-Analysis Architecture

Effect Size Estimation and the Winner's Curse

Following the discovery of significant rare variant associations, effect size estimation presents unique challenges due to both the winner's curse (upward bias of effect sizes in significant findings) and the heterogeneity of individual variant effects within aggregated sets [22]. The winner's curse arises because significant associations are more likely to be detected when the estimated effect size is larger than the true effect, particularly in underpowered settings [22]. For burden tests, this is further complicated by the fact that aggregated variant groups likely contain both causal and non-causal variants, with potentially bidirectional effects, leading to downward bias in the pooled effect estimate [22].

Several approaches have been proposed to address these biases, including bootstrap resampling methods and likelihood-based approaches that correct for the winner's curse [22]. However, these methods must contend with the complex interplay between upward bias from significance filtering and downward bias from effect heterogeneity. Research has shown that the magnitude and direction of bias depend on both the underlying genetic architecture and the specific association test used, with burden tests and variance-component tests exhibiting different bias patterns [22]. These considerations highlight the importance of interpreting reported effect sizes from rare variant studies with appropriate caution and using bias-correction methods when possible.

Table 3: Research Reagent Solutions for Rare Variant Studies

Research Tool	Function	Application Context
Exomiser	Variant prioritization and annotation	Identification of putative pathogenic variants from WES data
geneBurdenRD	R framework for burden testing	Rare disease gene association discovery in family-based designs
Meta-SAIGE	Rare variant meta-analysis	Scalable cross-cohort association testing
SAIGE-GENE+	Individual-level association testing	Gene-based tests accounting for relatedness and imbalance
SKAT/SKAT-O R package	Variance-component testing	Gene-based association with heterogeneous variant effects
PathVar	Pathogenic variant filtering pipeline	Bioinformatic prioritization of likely functional variants

The statistical rationale for gene-based burden testing of rare variants represents a paradigm shift from the single-variant approaches that have dominated genetic association studies for common variants. By aggregating signals across multiple rare variants within functional units, these methods successfully overcome the severe power limitations that would otherwise prevent the detection of rare variant contributions to complex traits and diseases [19] [20]. The continued development and refinement of burden tests, variance-component tests, and adaptive omnibus tests have provided researchers with a powerful toolkit for probing the role of rare variation in human phenotypes [21] [18].

Future methodological developments will likely focus on improving power through more sophisticated incorporation of functional annotations, enhancing methods for cross-ancestry rare variant analysis, and developing approaches that jointly model common and rare variants [20]. As biobank-scale sequencing resources continue to expand, with projects like the 100,000 Genomes Project and All of Us providing unprecedented sample sizes, rare variant association methods will play an increasingly central role in unraveling the genetic architecture of both rare and common diseases [18] [3]. The integration of these statistical advances with experimental functional validation will be essential for translating genetic discoveries into biological insights and therapeutic opportunities [10] [24].

A Practical Toolkit: Core Burden Tests and Advanced Methodologies

The exploration of genetic associations has expanded significantly beyond common variants to include rare variants, defined as those with a minor allele frequency (MAF) of less than 1-5% [25] [21]. While individually rare, these variants collectively contribute substantially to the genetic architecture of complex traits and diseases. However, their low frequency presents unique statistical challenges for association testing, as standard single-variant tests typically used in genome-wide association studies (GWAS) lack sufficient power due to the necessity for extremely large sample sizes to detect effects [25] [21]. This limitation has driven the development of region-based or gene-based association methods that aggregate information across multiple rare variants within a functional unit, such as a gene or pathway.

The methodological spectrum for rare variant analysis has evolved to include several distinct but complementary approaches. Burden tests, one of the earliest classes of methods, operate on the principle of collapsing multiple variants into a single genetic score [25] [17]. These tests assume that all rare variants in the target region influence the phenotype in the same direction and with similar effect magnitudes [25]. While powerful when this assumption holds true, burden tests suffer substantial power loss when many non-causal variants are included or when both risk-increasing and protective variants are present within the same region [25] [7].

Recognizing these limitations, variance-component tests like the Sequence Kernel Association Test (SKAT) were developed to allow for heterogeneous effect sizes and directions [25] [17]. SKAT builds upon the kernel machine regression framework and treats variant effects as random, testing whether their variance differs from zero [25]. This flexibility makes SKAT particularly powerful in the presence of protective and deleterious variants and null variants, though it becomes less powerful than burden tests when a large proportion of variants in a region are causal and act in the same direction [25].

To leverage the advantages of both approaches, hybrid tests such as SKAT-O were developed, which optimally combine burden and SKAT statistics within a unified framework [25] [7]. Simultaneously, p-value combination methods like the Aggregated Cauchy Association Test (ACAT) have emerged as powerful alternatives, especially when only a small number of variants within a region are truly causal [26] [7]. This methodological spectrum offers researchers a diverse toolkit for rare variant association discovery, with the choice of method dependent on the underlying genetic architecture of the trait, which is typically unknown a priori.

Statistical Foundations and Comparative Analysis

Methodological Principles

Burden Tests

Burden tests represent the foundational approach to rare variant association analysis. These methods operate by collapsing genetic information across multiple rare variants within a predefined region (typically a gene) into a single composite genotype score for each individual [17]. The most straightforward approach involves creating a binary indicator variable that scores an individual as 1 if they carry at least one rare allele across the variant set, and 0 otherwise [26]. More sophisticated approaches calculate a weighted sum of variants, with weights often dependent upon the alternate allele frequency (AAF) of variants [26]. This burden score is then tested for association with the phenotype of interest using regression frameworks.

The statistical model for burden tests typically assumes that each regression coefficient βj in the generalized linear model can be expressed as βj = wjβ0, where wj is a predetermined weight for variant j (often based on MAF), and β0 is a common effect size [25]. This reduces the multiple testing problem from p degrees of freedom (where p is the number of variants) to a single degree of freedom test for H0: β0 = 0 [25]. While this collapsing strategy increases power when all variants are causal and have effects in the same direction, it suffers from notable limitations when these assumptions are violated, particularly when non-causal variants are included or when effects are bidirectional [25] [7].

Sequence Kernel Association Test (SKAT)

SKAT takes a fundamentally different approach by treating the variant effects as random rather than fixed. Specifically, it assumes that each βj independently follows an arbitrary distribution with mean zero and variance ψwj², where wj is a fixed weight for variant j [25]. Under this assumption, testing the null hypothesis H0: β = 0 is equivalent to testing H0: ψ = 0, which represents a variance component test in generalized linear mixed models [25].

The SKAT test statistic is derived from the score statistic of this variance component and takes the form of a quadratic function of the residuals from the null model (without genetic effects) [25]. For continuous traits, this can be expressed as:

[ Q = (y - \hat{\mu})'K(y - \hat{\mu}) ]

where K = GW²G' is an n×n weighted linear kernel matrix, G is the n×p genotype matrix, W is a p×p diagonal matrix of weights, and (y - μ̂) is the residual vector from the null model [25]. SKAT's key advantage is its flexibility—it does not assume all variants have effects in the same direction, making it robust to the presence of both risk and protective variants in the same gene [25] [17].

SKAT-O: The Optimal Unified Test

SKAT-O was developed to bridge the gap between burden tests and SKAT by creating an optimal test that adapts to the underlying genetic architecture [25]. This method considers a class of tests that are linear combinations of the burden test and SKAT statistics, and identifies the test that maximizes power [25]. Mathematically, SKAT-O can be formulated as a generalized family of SKAT tests that incorporates a correlation structure of variant effects through a family of kernels [25].

The method estimates a correlation parameter ρ in the kernel matrix that determines the balance between burden and SKAT components, with ρ = 1 corresponding to a burden test and ρ = 0 corresponding to SKAT [25]. The test statistic is based on a mixture distribution that depends on this parameter, and p-values can be calculated analytically without requiring computationally intensive resampling procedures, making it feasible for genome-wide applications [25]. Simulation studies have demonstrated that SKAT-O outperforms both burden tests and SKAT across a wide range of scenarios, particularly when the proportion of causal variants and their effect directions are unknown [25].

Aggregated Cauchy Association Test (ACAT)

ACAT represents a different philosophical approach to rare variant testing by combining p-values from single-variant tests or other complementary association tests [26] [7]. This method transforms individual p-values into Cauchy variables and aggregates them into a single test statistic [7]. The key advantage of ACAT is its robustness and power when only a small number of variants in a region are truly causal [26] [7].

ACAT is particularly useful for combining results across different functional annotations or allele frequency bins, which has become common practice in comprehensive rare variant analyses [26] [27]. For example, researchers often test multiple variant sets for each gene, such as loss-of-function (LoF) variants only, LoF plus missense variants, and different MAF cutoffs (≤1%, ≤0.1%, ≤0.01%) [27]. ACAT can efficiently combine these results while maintaining appropriate type I error control and computational efficiency [26].

Comparative Performance Characteristics

Table 1: Comparative analysis of rare variant association tests

Method	Key Assumptions	Strengths	Limitations	Optimal Use Case
Burden Tests	All variants are causal with same direction and similar effect sizes [25]	High power when assumptions hold; simple interpretation [25] [21]	Power loss with non-causal variants or bidirectional effects [25] [7]	Genes with clear LoF mechanism; all variants deleterious [21]
SKAT	Variant effects random with mean zero; allows bidirectional effects [25]	Robust to mixed effect directions; no need to specify causal variants [25] [17]	Less power than burden when high proportion of causal variants with same direction [25]	Genomic regions with both risk and protective variants [25]
SKAT-O	Adaptive method that bridges burden and SKAT [25]	Optimal across wider range of scenarios; data-adaptive [25]	Computational complexity; requires estimating correlation parameter [25]	Default choice when genetic architecture unknown [25]
ACAT	Combines evidence from multiple tests [26] [7]	Powerful with few causal variants; combines different annotation/frequency bins [26] [7]	May lose power with many non-causal variants [7]	Combining multiple masks; small number of causal variants [26] [7]
Sparteine	Sparteine|CAS 90-39-1|Research Chemical	Sparteine is a quinolizidine alkaloid for research use only (RUO). Explore its applications in neuroscience, antiarrhythmic studies, and as a chiral ligand. Not for human consumption.	Bench Chemicals
(-)-Vasicine	(-)-Vasicine, CAS:6159-55-3, MF:C11H12N2O, MW:188.23 g/mol	Chemical Reagent	Bench Chemicals

The performance of these methods depends critically on the underlying genetic architecture, particularly the proportion of causal variants and the consistency of their effect directions. Empirical and analytical studies have shown that burden tests are more powerful than single-variant tests only when a substantial proportion of variants are causal, while variance component tests like SKAT maintain power across more diverse scenarios [21].

Recent methodological developments have focused on enhancing these core approaches. The Sparse Burden Association Test (SBAT) jointly models a set of burden scores under the assumption that causal burden scores act in the same effect direction, simultaneously assessing significance and selecting the burden scores that best explain the association [26]. SAIGE-GENE+ addresses computational and accuracy limitations in large-scale biobank data by collapsing ultra-rare variants (minor allele count ≤ 10) before testing, greatly improving type I error control for very rare variants while maintaining power [27].

Application Notes and Experimental Protocols

Protocol 1: Gene-Based Rare Variant Association Analysis

Study Design and Sample Size Considerations

An effective rare variant association study requires careful consideration of sample size, which depends on region heritability (h²), the number of causal variants (c), and the total number of variants tested (v) [21]. Analytical calculations and simulations indicate that aggregation tests are more powerful than single-variant tests only when a substantial proportion of variants are causal [21]. For example, when aggregating all rare protein-truncating variants (PTVs) and deleterious missense variants, aggregation tests become more powerful than single-variant tests for >55% of genes when PTVs, deleterious missense variants, and other missense variants have 80%, 50%, and 1% probabilities of being causal, with n=100,000 and h²=0.1% [21].

Variant Annotation and Quality Control

The first step involves comprehensive variant annotation and quality control. For whole exome or genome sequencing data, implement the following QC pipeline:

• Variant Calling and Filtering: Begin with raw sequencing data (FASTQ files) processed through a standardized variant calling pipeline (e.g., GATK best practices). Filter variants based on call rate (>95%), Hardy-Weinberg equilibrium (p > 1×10⁻⁶), and sequencing depth.

• Functional Annotation: Annotate variants using tools like ANNOVAR or VEP (Variant Effect Predictor) to classify them by functional consequence (synonymous, missense, LoF, etc.). Prioritize potentially functional variants including protein-truncating variants (PTVs), predicted deleterious missense variants (e.g., REVEL score > 0.5), and splice region variants [26] [27].

• Population Frequency Filtering: Determine variant frequencies using large reference populations (gnomAD, UK Biobank). Define rare variants based on MAF threshold (typically < 0.5-1%) [21]. Consider multiple frequency bins (MAF ≤ 1%, ≤ 0.1%, ≤ 0.01%) to capture frequency-dependent effects [27].

• Variant Set Definition: Define variant sets for testing, typically at the gene level, but consider functional sub-units or domains for large genes. Create multiple annotation categories: (1) LoF only; (2) LoF + deleterious missense; (3) LoF + all missense [27].

Association Testing Workflow

The following workflow implements a comprehensive rare variant association analysis:

• Null Model Fitting: For each phenotype, fit a null model including relevant covariates (age, sex, principal components to account for population stratification). For binary traits with unbalanced case-control ratios, use methods that account for this imbalance, such as SAIGE-GENE+ [27].

• Single-Variant Testing: Perform single-variant association tests as a baseline. While underpowered for individual rare variants, these provide complementary information and can be combined using ACAT-V [26].

• Gene-Based Testing: Conduct the spectrum of gene-based tests:

  • Burden tests using weighted sums of variants (weights often based on MAF using beta density function) [25]
  • SKAT using linear kernels [25]
  • SKAT-O to optimally combine burden and SKAT approaches [25]
  • ACAT to combine p-values across different functional annotations and frequency bins [26] [27]

• Multiple Testing Correction: Apply gene-based multiple testing correction, with exome-wide significance threshold typically set at p < 2.5×10⁻⁶ (0.05/20,000 genes) [27].

• Validation and Replication: Seek replication in independent cohorts where possible. For novel findings, conduct functional validation experiments to establish biological mechanisms.

The following diagram illustrates the key decision points in selecting appropriate statistical tests based on the expected genetic architecture:

Protocol 2: Application to Large-Scale Biobank Data

The emergence of large-scale biobanks with whole exome or genome sequencing data requires specialized analytical approaches to address computational challenges and specific characteristics of biobank data.

SAIGE-GENE+ Implementation for Biobank Data

SAIGE-GENE+ was specifically developed to address limitations in analyzing biobank-scale sequencing data [27]. Implement the following protocol:

• Data Preparation: Organize genotype data in efficient formats (BGEN, PLINK). Prepare phenotype and covariate files, ensuring consistent sample identifiers.

• Step 1: Null Model Fitting: Fit the null logistic mixed model for binary traits or linear mixed model for quantitative traits, including a genetic relationship matrix (GRM) to account for population structure and relatedness:

• Step 2: Association Testing: Perform set-based rare variant tests, specifying multiple MAF cutoffs and functional annotations:

• Ultra-Rare Variant Handling: SAIGE-GENE+ automatically collapses variants with minor allele count (MAC) ≤ 10 before testing to address data sparsity, which improves type I error control for very rare variants [27].

• Result Integration: Use ACAT to combine results across different MAF cutoffs and functional annotations to maximize power [27].

Sparse Burden Association Test (SBAT) for Multiple Burden Scores

When testing multiple correlated burden scores across different annotation classes and frequency thresholds, SBAT provides an effective approach:

• Burden Score Calculation: Compute multiple burden scores for each gene across different AAF thresholds (e.g., ≤1%, ≤0.1%, ≤0.01%) and annotation classes (pLoFs, deleterious missense) [26].

• Joint Model Fitting: Fit a joint model of all burden scores using non-negative least squares (NNLS) with positivity constraints on regression parameters, enforcing the prior distribution that causal burden scores act in the same direction [26].

• Statistical Testing: Calculate the quadratic form test statistic T = β̂'V⁻¹β̂, where β̂ is the NNLS estimate and V is the covariance matrix. The null distribution follows a mixture of chi-square distributions [26].

• Two-Tailed Testing: Apply SBAT twice to both Y and -Y to account for unknown effect direction, then combine p-values using the Cauchy combination method [26].

Table 2: Key computational tools and resources for rare variant association studies

Tool/Resource	Primary Function	Application Notes	Reference
REGENIE	Whole-genome regression for polygenic trait analysis	Implements burden tests, SKAT, SKAT-O, ACAT; efficient for large datasets	[26]
SAIGE-GENE+	Set-based rare variant tests for biobank data	Addresses case-control imbalance; collapses ultra-rare variants (MAC ≤ 10)	[27]
geneBurdenRD	R framework for rare disease gene burden testing	Specialized for Mendelian diseases; handles unbalanced case-control studies	[3]
Exomiser	Variant prioritization tool	Filters and annotates putative causal variants; integrates phenotypic data	[3]
PLINK/SEQ	Toolset for genome association analysis	Implements various rare variant association tests; flexible variant filtering	[7]
ANNOVAR/VEP	Functional variant annotation	Annotates consequences, population frequencies, pathogenicity predictions	[3]
UK Biobank WES	Reference dataset with exome sequencing	160,000 individuals; 30 quantitative/141 binary traits for calibration	[27]
gnomAD	Population frequency reference	Annotates variant frequencies across diverse populations	[21]

The methodological spectrum for rare variant association analysis has evolved substantially, offering researchers a diverse toolkit tailored to different genetic architectures and study designs. Burden tests remain powerful when the underlying biology suggests unidirectional effects across variants, while SKAT provides robustness to heterogeneous effect directions. SKAT-O offers an adaptive approach that maintains power across diverse scenarios, and ACAT excels when combining evidence across multiple tests or when only a few variants are causal.

Future methodological developments will likely focus on integrating multi-omics data, improving cross-population portability, and enhancing statistical power for very rare variants through sophisticated collapsing methods and functional priors. As sequencing studies continue to grow in scale, efficient implementations like SAIGE-GENE+ and REGENIE will become increasingly important for analyzing biobank-scale data while maintaining appropriate type I error control. The ongoing challenge remains the unknown genetic architecture for most traits, necessishing continued development of robust, adaptive methods and thoughtful application of the existing methodological spectrum.

Gene-based burden tests have become a cornerstone of rare variant association studies in rare disease research and drug target discovery. These tests operate on a fundamental principle: aggregating the effects of multiple rare genetic variants within a gene to increase statistical power for detecting associations with diseases or traits. The core of this methodology hinges on the accurate selection and annotation of variant types deemed most likely to be functional—primarily protein-truncating variants (pLoF) and deleterious missense variants. The strategic definition of these categories and the implementation of appropriate frequency thresholds directly determine the success of gene-based burden tests in distinguishing true biological signals from background noise. This protocol details the standardized frameworks and computational methodologies required for rigorous variant selection and annotation, specifically tailored for gene-based burden analyses in rare variant research.

Defining Key Variant Classes

Protein-Truncating Variants (pLoF)

Protein-truncating variants (pLoFs) represent a category of high-impact genetic alterations that prematurely terminate protein translation, typically resulting in a loss of function of the affected gene. pLoF variants are considered the most severe in their predicted functional consequences and are often prioritized in rare variant burden analyses due to their high prior probability of pathogenicity. pLoF variants include nonsense variants (single nucleotide substitutions that create a premature stop codon), frameshift variants (insertions or deletions that disrupt the reading frame, leading to premature termination), and canonical splice-site variants (mutations at highly conserved intronic positions that disrupt mRNA splicing, often resulting in nonsense-mediated decay). Research indicates pLoF variants have a high probability (approximately 80%) of being causal in disease contexts [21].

Missense Variants

Missense variants are single nucleotide substitutions that result in the substitution of one amino acid for another in the translated protein. The functional impact of missense variants exists on a continuous spectrum, ranging from completely neutral to severely deleterious. This makes their interpretation particularly challenging. In burden testing, the critical task is to distinguish deleterious missense variants (those likely to impair protein function and contribute to disease) from neutral variants. This discrimination relies heavily on computational prediction algorithms and functional assessments. Current evidence suggests that putatively deleterious missense variants have an estimated 50% probability of being causal, while other missense variants have a much lower probability (approximately 1%) [21]. This distinction is crucial for constructing powerful burden test masks.

Frequency Thresholds for Rare Variants

Frequency thresholds are applied to filter variants based on their population allele frequency, focusing the analysis on rare variation. The standard threshold for defining a rare variant is a Minor Allele Frequency (MAF) < 1% [21] [18]. However, more stringent thresholds (e.g., MAF < 0.1% or MAF < 0.01%) are often employed for gene-based burden tests to further enrich for very recent or severe mutations. The selection of a specific threshold involves a trade-off between signal purity and the number of variants available for analysis. For extremely rare variants, some methods like Meta-SAIGE collapse ultrarare variants (MAC < 10) to boost power and improve type I error control [18].

Table 1: Key Variant Classes and Their Characteristics in Burden Testing

Variant Class	Molecular Consequence	Typical Prioritization in Burden Tests	Estimated Probability of Being Causal
Protein-Truncating (pLoF)	Premature termination of protein translation	Highest priority	~80% [21]
Deleterious Missense	Amino acid change predicted to impair function	High priority	~50% [21]
Other Missense	Amino acid change with unknown/neutral impact	Often excluded or down-weighted	~1% [21]
Synonymous	No change to amino acid sequence	Typically excluded	Very low

Computational Annotation and Prediction Methods

Protein Language Models for Missense Variant Effect Prediction

Protein Language Models (PLMs), such as ESM1b and ESM2, are deep learning models pre-trained on millions of protein sequences, allowing them to learn the underlying principles of protein structure and function. These models can predict the functional impact of missense variants without relying on multiple sequence alignments, providing a fast and scalable annotation method. The ESM1b model calculates a log-likelihood ratio (LLR) score, representing the model's assessment of the plausibility of the mutant amino acid in its structural context compared to the wild-type. More negative scores indicate a higher probability of the variant being deleterious. ESM1b has demonstrated superior performance in classifying pathogenic and benign variants in benchmarks like ClinVar and HGMD, outperforming many existing methods [28]. PLMs can be further fine-tuned for specific tasks, such as predicting functional protein features (e.g., active sites, binding residues) at amino-acid resolution, adding a layer of mechanistic interpretation to variant impact [29].

Integrated Annotation Pipelines and Tools

For practical implementation, researchers often use integrated bioinformatics pipelines that consolidate multiple annotation sources. Tools like Exomiser and Genomiser perform integrated variant prioritization by combining genotype, pathogenicity prediction, allele frequency, and phenotype data (encoded using Human Phenotype Ontology/HPO terms) to generate a ranked list of candidate variants or genes [15]. The PathVar pipeline is an example of a tool designed specifically to identify pathogenic variants for inclusion in burden tests, and its use has been shown to contribute uniquely to identifying gene associations in complex disorders like hemiplegic migraine [4]. These pipelines typically incorporate a meta-prediction score from dbNSFP (Database for Nonsynonymous SNPs' Functional Predictions), which aggregates results from dozens of functional prediction algorithms [28].

Diagram 1: Workflow for variant selection and annotation prior to gene-based burden testing. Key decision points include frequency filtering and functional annotation.

Protocol for Variant Selection for Gene-Based Burden Testing

Variant Quality Control and Frequency Filtering

Begin with high-quality variant calls in VCF format, derived from exome or genome sequencing. Apply standard quality control (QC) filters, such as call rate (>95%), genotype quality, and deviation from Hardy-Weinberg equilibrium. Subsequently, filter variants based on population frequency using large, diverse reference databases like gnomAD. The recommended threshold for a standard rare variant analysis is MAF < 1%. For studies aiming to capture ultra-rare, high-effect variants, a more stringent threshold of MAF < 0.1% or 0.01% is advised. This step is critical for excluding common polymorphisms unlikely to cause severe rare diseases.

Functional Annotation and Mask Definition

Annotate the QC-passed rare variants using the following steps:

• pLoF Annotation: Use tools like LOFTEE or VEP to confidently identify high-confidence pLoF variants, excluding those in low-confidence regions or near the end of transcripts.
• Missense Impact Prediction: Score all missense variants using a high-performance predictor such as ESM1b. The recommended threshold for classifying a variant as "deleterious" is an ESM1b LLR score < -7.5, which corresponds to a true-positive rate of 81% and a true-negative rate of 82% for known pathogenic variants [28].
• Construct a Mask: Define the set (or "mask") of variants to be included in the burden test. A common and effective strategy is to create a mask that aggregates all pLoF variants and deleterious missense variants (defined by the ESM1b threshold) within a gene [21].

Gene-Based Burden Test Implementation

With the variant mask defined, proceed with the gene-based burden test. For case-control studies, tools like SAIGE-GENE+ are recommended as they effectively handle case-control imbalance and sample relatedness [18]. In scenarios where only case data is available (case-only designs), methods like the Case-Only Burden Test (COBT) can be applied. COBT uses a Poisson model to test for an excess of damaging variants in a gene compared to expectations based on general population mutation rates [5]. For meta-analysis across multiple cohorts, Meta-SAIGE provides a scalable solution that controls type I error effectively and achieves power comparable to pooled analysis of individual-level data [18].

Table 2: Key Tools and Resources for Variant Selection and Burden Testing

Tool/Resource	Primary Function	Application in Protocol	URL/Reference
gnomAD	Population frequency database	Frequency filtering (MAF < 1%)	gnomad.broadinstitute.org
ESM1b/ESM2	Protein Language Model	Missense variant effect prediction (LLR score)	[28] [29]
Exomiser	Integrated variant prioritization	Combines genotype, phenotype & prediction for ranking	github.com/exomiser/Exomiser [15]
SAIGE-GENE+	Rare variant association test	Gene-based burden testing in case-control cohorts	[18]
COBT	Case-Only Burden Test	Gene-based burden testing in case-only cohorts	github.com/RausellLab/COBT [5]
Meta-SAIGE	Rare variant meta-analysis	Combining burden test results across cohorts	[18]

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Computational Tools and Databases for Variant Burden Studies

Item	Specifications/Version	Function in Experiment
ESM1b Model	650 million parameters	Provides unsupervised missense variant effect scores (LLR) for classifying deleterious variants. Outperforms many existing methods [28].
Exomiser Suite	v13.2.0+	Open-source software for phenotype-driven variant prioritization. Combines frequency, pathogenicity predictions, and HPO-based gene-phenotype associations [15].
ClinVar Database	Latest public release	Public archive of reported variant-disease relationships. Serves as a key benchmark for validating pathogenicity predictions [28] [29].
GRCh38/hg38	Human reference genome	The current standard reference genome for alignment and variant calling, essential for accurate annotation.
SAIGE-GENE+ Software	v1.2+	Statistical tool for performing gene-based rare variant tests (Burden, SKAT, SKAT-O) while accounting for sample relatedness and imbalance [18].
1000 Genomes Project	Phase 3 data	Serves as a source of control genotypes and for estimating population-specific allele frequencies.
TGX-155	TGX-155, CAS:351071-90-4, MF:C20H19FN2O3, MW:354.4 g/mol	Chemical Reagent
Wilforlide A	Wilforlide A, CAS:84104-71-2, MF:C30H46O3, MW:454.7 g/mol	Chemical Reagent

Advanced Considerations and Troubleshooting

Optimizing Parameters for Diagnostic Prioritization

When using tools like Exomiser for diagnostic variant prioritization, parameter optimization is critical. Evidence from analyses of the Undiagnosed Diseases Network (UDN) suggests that optimizing parameters—such as the choice of gene-phenotype association data and variant pathogenicity predictors—can significantly improve performance. For instance, optimized parameters increased the ranking of coding diagnostic variants within the top 10 candidates from 67.3% to 88.2% in exome sequencing data [15]. It is also important to consider isoform-specific variant effects, as a variant may be damaging in one protein isoform but neutral in another. Protein language models like ESM1b can predict these isoform-sensitive effects, improving annotation accuracy [28].

When are Burden Tests Most Powerful?

Burden tests are not universally the most powerful approach. Their performance relative to single-variant tests depends heavily on the underlying genetic architecture. Analytical calculations and simulations reveal that burden tests are more powerful than single-variant tests only when a substantial proportion of the aggregated variants are causal and have effects in the same direction [21]. For example, if pLoF and deleterious missense variants have high probabilities of being causal (e.g., 80% and 50%, respectively), burden tests are more powerful for >55% of genes. However, if the set of aggregated variants is diluted with many neutral variants, single-variant tests may be superior.

Diagram 2: Decision framework for choosing between burden tests and single-variant tests based on genetic architecture.

Addressing Interpretation Challenges for VUS

A significant challenge in clinical interpretation is the reclassification of Variants of Uncertain Significance (VUS). Protein language models fine-tuned on specific functional features (e.g., active sites, binding domains) can provide mechanistic insights by predicting if a VUS disrupts a critical protein feature [29]. Furthermore, experimental data from multiplexed assays of variant effect (MAVEs), which systematically measure the functional impact of thousands of variants in a single experiment, can be used as a gold standard for clinical interpretation. For instance, a MAVE of SOD1 provided functional evidence for reclassifying 41% of VUSs associated with amyotrophic lateral sclerosis [30]. Integrating these functional maps with computational predictions is a powerful future direction for variant annotation.

The pursuit of understanding the genetic underpinnings of complex diseases and traits has increasingly shifted focus towards the role of rare genetic variants. These variants, typically defined as having a Minor Allele Frequency (MAF) of less than 1%, are thought to have substantial effects on disease risk and phenotypic expression. However, their individual rarity presents a significant statistical challenge, as traditional single-variant association tests lack the power to detect their effects reliably. To overcome this, gene-based burden tests have been developed as a powerful strategy. These methods aggregate, or "collapse," multiple rare variants within a predefined genomic unit—such as a gene—into a single composite test statistic, thereby amplifying the signal for association. Within the broader thesis of rare variant research, this approach allows investigators to test the collective contribution of multiple functionally related variants to a phenotype, providing a more holistic view of gene function and dysfunction.

Several sophisticated statistical software tools have been engineered to perform these tests efficiently on large-scale genomic data. This protocol details the application of three prominent methods: the TRAPD (Trans-Omics for Precision Medicine Whole Genome Sequencing Data Analysis Plan) design principle, and the software implementations SAIGE-GENE+ and REGENIE. The TRAPD framework provides a robust standard for the meta-analysis of rare variant associations across multiple cohorts. SAIGE-GENE+ and REGENIE are cutting-edge tools designed to handle the computational and statistical complexities of burden testing in biobank-scale datasets, particularly for case-control phenotypes with extreme imbalance. This article provides detailed Application Notes and Protocols for researchers, scientists, and drug development professionals seeking to implement these tests, complete with structured data comparisons, experimental workflows, and essential reagent solutions.

Table 1: Key Characteristics of Featured Burden Testing Tools

Tool/Method	Primary Use Case	Core Statistical Approach	Handling of Case-Control Imbalance	Relatedness Adjustment
TRAPD Design	Multi-cohort meta-analysis	Aggregation of summary statistics	Addressed by cohort-level analysis	Uses genetic relationship matrix (GRM) or polygenic effects [31]
SAIGE-GENE+	Set-based rare variant tests in large biobanks	Burden, SKAT, SKAT-O	Saddle Point Approximation (SPA)	Full or sparse GRM [32]
REGENIE	Whole-genome regression for common & rare variants	Burden tests, SKAT, SKAT-O, ACAT	Firth logistic regression, SPA	Polygenic effect estimates [33]

Tool Comparison and Selection Guidelines

Selecting the appropriate tool for a rare variant burden analysis is critical and depends on factors such as sample size, data structure, and the specific biological question. The following section provides a comparative overview of the featured tools to guide this decision.

SAIGE-GENE+ is a specialized tool for set-based rare variant association tests that excels in large datasets with highly unbalanced case-control ratios, a common scenario in biobank studies of rare diseases. Its key innovation is the use of Saddle Point Approximation (SPA) to accurately calibrate p-values in the presence of extreme imbalance, which prevents inflation of false positive rates. A significant improvement over its predecessor, SAIGE-GENE+ collapses ultra-rare variants (with Minor Allele Count, MAC ≤ 10) to mitigate data sparsity, leading to superior type I error control and a dramatic reduction in computation time and memory usage. For instance, testing the large gene TTN required 1,407 times less CPU time and reduced memory usage from 65 GB to 2.1 GB compared to SAIGE-GENE [32]. SAIGE-GENE+ facilitates powerful analyses by allowing multiple tests per gene using different MAF cutoffs and functional annotations.

REGENIE offers a unified framework for both common and rare variant association studies. It employs a two-step whole-genome regression approach that is highly computationally efficient for massive datasets. For burden testing, it can implement various statistical models, including Burden, SKAT, and ACAT. REGENIE controls for population and family structure using polygenic effect estimates rather than a direct GRM. However, independent evaluations suggest that REGENIE might not be the optimal choice for analyzing highly correlated data of small size, such as large pedigrees, where other tools like GENESIS or SAIGE might be more robust [31].

The TRAPD framework is not a single software tool but a design and analysis standard for multi-cohort rare variant association studies. It outlines a protocol for performing cohort-level burden tests using robust methods (which could include SAIGE-GENE+ or REGENIE) followed by a meta-analysis that aggregates summary statistics across studies. This approach is powerful for increasing sample size and discovering associations that are reproducible across diverse populations.

Table 2: Performance and Practical Considerations for Burden Testing Tools

Aspect	SAIGE-GENE+	REGENIE	Notes
Computational Efficiency	High (after improvements)	High	SAIGE-GENE+ memory usage dropped to 2.1 GB for large genes [32].
Optimal Data Structure	Large biobanks, unbalanced case-control	Large datasets, less correlated data	REGENIE may underperform in small, highly correlated datasets [31].
Type I Error Control	Excellent with SPA and collapsing	Good, but can be conservative	fastGWA-GLMM, another tool, was noted as overly conservative in family data [33].
Key Advantage	Superior control for extreme imbalance; integrated multiple testing	Unified GWAS and rare variant pipeline; speed	SAIGE-GENE+ identified 551 gene-phenotype associations in UK Biobank data [34].

Experimental Protocol for Burden Testing

This protocol outlines a standardized workflow for conducting gene-based burden tests, from data preparation to result interpretation. The steps can be adapted for use with SAIGE-GENE+, REGENIE, or within the TRAPD meta-analysis framework.

Pre-processing and Quality Control of Genetic Data

• Variant Calling and Initial QC: Begin with high-quality whole genome or whole exome sequencing data that has been aligned, variant-called (e.g., using GATK), and subjected to initial quality control. This includes removing samples with excessive missingness, checking for sex discrepancies, and eliminating variants that are monomorphic or have a low call rate (e.g., < 90%) [31].
• Variant Annotation: Annotate the remaining variants using a tool like Ensembl VEP or SnpEff to identify their functional consequences (e.g., LoF, missense, synonymous).
• Defining Qualifying Variants (QVs): For each gene, define the set of rare variants to be aggregated. A common strategy is to include:
  • LoF variants: High-confidence loss-of-function variants (e.g., stop-gain, stop-loss, splice-site).
  • Damaging missense variants: Missense variants predicted to be deleterious by tools like PolyPhen-2 or SIFT.
  • MAF Filtering: Restrict variants based on a MAF threshold (e.g., < 0.1%, < 1%) within the study population. Using multiple MAF cutoffs can enhance power, as associations are sometimes enriched in very rare variants [32].
• Sample and Covariate File Preparation: Prepare a phenotype file containing the case-control status or quantitative trait values. Prepare a covariate file containing relevant adjustments such as age, sex, sequencing batch, and genetic principal components (PCs) to account for population stratification.

Running Association Tests with SAIGE-GENE+

SAIGE-GENE+ operates in two steps [32]:

Step 1: Fitting the Null Model The tool first fits a null logistic mixed model containing only the covariates and a random effect to account for relatedness. This step calculates the genetic relationship matrix (GRM) from all available genotypes.

Step 2: Performing Association Tests In the second step, the model residuals from Step 1 are used to test for association between the aggregated variant sets and the phenotype. SAIGE-GENE+ can perform Burden, SKAT, and SKAT-O tests simultaneously.

Running Association Tests with REGENIE

REGENIE also uses a two-step approach [35] [33]:

Step 1: Fitting the Null Model REGENIE first fits a null model using a whole-genome regression technique to estimate polygenic effects, which capture population structure and relatedness.

Step 2: Association Testing In Step 2, REGENIE tests for associations, conditioning on the polygenic effects from Step 1. For burden tests, variants are collapsed within each gene.

Meta-Analysis Following the TRAPD Framework

For multi-cohort studies, the TRAPD design involves:

• Cohort-Level Analysis: Each participating cohort runs gene-based burden tests locally using a standardized protocol (e.g., using SAIGE-GENE+ with identical QC steps, MAF cutoffs, and variant annotations).
• Summary Statistics Sharing: Each cohort prepares a file of summary statistics for each gene test, including the effect size (beta), standard error, p-value, and other relevant metrics.
• Meta-Analysis: A central coordinating center performs a fixed- or random-effects meta-analysis (using tools like METAL or R's metafor package) to combine the summary statistics across all cohorts, generating a single, combined association estimate for each gene-phenotype pair.

Visualization of the Burden Testing Workflow

The following diagram, generated using Graphviz, illustrates the logical flow and decision points in a standard gene-based burden testing protocol, integrating the roles of SAIGE-GENE+, REGENIE, and the TRAPD framework.

Diagram 1: A generalized workflow for gene-based burden testing, showing the key steps from data preparation through analysis with SAIGE-GENE+ or REGENIE, and integration into a TRAPD meta-analysis for multi-cohort studies.

The Scientist's Toolkit: Research Reagent Solutions

The following table details the essential software and data "reagents" required to implement the burden testing protocols described in this article.

Table 3: Essential Research Reagents for Gene-Based Burden Tests

Reagent / Resource	Type	Function / Application	Example or Note
SAIGE-GENE+	Software Tool	Performs set-based rare variant association tests (Burden, SKAT, SKAT-O) for binary traits with imbalance [32].	https://github.com/weizhouUMICH/SAIGE
REGENIE	Software Tool	Performs whole-genome regression for common and rare variant association tests, including burden tests [35].	https://github.com/rgcgithub/regenie
GATK	Software Tool	Industry standard for variant discovery in high-throughput sequencing data; used for initial variant calling [31].	Best practices pipeline includes HaplotypeCaller.
ENSEMBL VEP	Software Tool	Annotates variants with functional consequences (e.g., LoF, missense), critical for defining qualifying variants [32].	Alternative: SnpEff.
Genetic Relationship Matrix (GRM)	Data / Model	A matrix quantifying genetic similarity between individuals, used to control for relatedness and population structure [31].	Calculated from genome-wide genotypes; used in SAIGE-GENE+ null model.
Phenotype & Covariate File	Data File	Tab-delimited file containing the trait of interest and adjustment variables (e.g., age, sex, principal components).	Essential input for all association testing tools.
Gene Annotation File	Data File	File (e.g., GTF format) defining the genomic coordinates of genes, used to group variants for burden tests.
UK Biobank WES Data	Reference Data	Large-scale, publicly available exome sequencing data used for method development and application [32].	SAIGE-GENE+ identified 551 gene-phenotype associations in this data.
Toddalolactone	Toddalolactone, CAS:483-90-9, MF:C16H20O6, MW:308.33 g/mol	Chemical Reagent	Bench Chemicals
Squalene	Squalene, CAS:111-02-4, MF:C30H50, MW:410.7 g/mol	Chemical Reagent	Bench Chemicals

Gene-based burden tests have emerged as a powerful statistical framework for identifying associations between rare genetic variants and complex phenotypes. By aggregating rare variants within a gene, these methods amplify the signal from individual variants that would otherwise be too rare to detect in single-variant association analyses [21]. This approach has proven particularly valuable in cardiometabolic research, where it has uncovered novel genetic determinants of diseases and biomarkers, providing insights into biological mechanisms and potential therapeutic targets. This application note presents key success stories and detailed methodologies from the application of gene-based burden testing to cardiometabolic traits and blood biomarkers, framing them within the broader context of rare variant research.

Success Stories in Gene-Based Burden Testing

Novel Gene Discoveries for Cardiometabolic Diseases

Large-scale genomic studies applying gene-based burden tests have successfully identified several previously unknown genetic associations with cardiometabolic diseases. These discoveries highlight the power of rare variant aggregation in uncovering biologically significant pathways.

Table 1: Novel Gene Associations Discovered Through Burden Testing

Gene	Phenotype	Variant Type	Effect Size	Biological Mechanism
IRS2	Type 2 Diabetes	Protein-truncating variants	OR = 6.4 [36]	Insulin signaling pathway component
RIF1	Body Mass Index	Protein-truncating variants	β = 2.66 kg/m² [36]	DNA double-strand break repair
UBR3	BMI & T2D	Protein-truncating variants	Increased risk for both traits [36]	E3-ubiquitin ligase activity
LDLR	LDL Cholesterol	Rare deleterious variants	Positive effect (z-score: 23.4) [17]	LDL receptor function
PCSK9	LDL Cholesterol	Rare deleterious variants	Negative effect (z-score: -14.2) [17]	Regulation of LDL receptor degradation

The discovery of protein-truncating variants in IRS2 (Insulin Receptor Substrate 2) represents a particularly significant finding. IRS2 encodes a key adapter molecule in the insulin signaling cascade, and rare damaging variants in this gene were found to confer a 6.4-fold increased risk for type 2 diabetes, with 34% of carriers developing the disease [36]. This finding not only validates IRS2 as a critical player in human glucose homeostasis but also suggests that impairment of IRS2-mediated signaling may represent an important mechanism in diabetes pathogenesis.

Similarly, the identification of RIF1 variants associated with BMI highlights how burden testing can reveal unexpected biological connections. RIF1 primarily functions in DNA double-strand break repair through the non-homologous end-joining pathway, yet protein-truncating variants were associated with substantially higher BMI (β = 2.66 kg/m²) [36]. This suggests potential novel connections between DNA repair mechanisms and metabolic regulation that warrant further investigation.

Blood Biomarker Regulation Through Rare Variants

Gene-based burden tests have also elucidated the genetic architecture of quantitative blood biomarkers, revealing how rare variants contribute to their regulation.

Table 2: Blood Biomarker Associations with Rare Variants

Biomarker	Significantly Associated Genes	Notable Example	Effect Direction
Alkaline Phosphatase	ALPL, others	ALPL	Negative (z-score: -49.6) [17]
LDL Direct Measurement	LDLR, PCSK9, NPC1L1, ABCG8	LDLR	Positive (z-score: 23.4) [17]
Sex Hormone Binding Globulin	SHBG	SHBG	Negative [17]
Testosterone	SHBG	SHBG	Negative [17]

The association between ALPL and alkaline phosphatase levels demonstrates how burden tests can identify genes with extremely strong effects on biomarkers. Rare functional variants in ALPL showed the largest effect size (z-score: -49.6) across all biomarker analyses, reflecting that ALPL encodes the alkaline phosphatase protein itself, and damaging variants directly reduce its production and activity [17].

The contrasting effects on LDL cholesterol regulation illustrate how burden tests can capture complementary biological pathways. While rare deleterious variants in LDLR (LDL receptor) increase LDL levels due to reduced clearance, variants in PCSK9 show the opposite effect, with damaging variants resulting in lower LDL levels due to increased LDL receptor availability [17]. This latter effect mirrors the mechanism of PCSK9 inhibitor drugs, validating the approach for identifying potential therapeutic targets.

Experimental Protocols

Gene-Based Burden Test Workflow

The standard analytical workflow for gene-based burden tests involves multiple stages from variant selection to statistical testing and validation. The following diagram illustrates this process:

Gene Burden Test Workflow

Variant Quality Control and Annotation

Begin with raw sequencing data (whole genome or exome) and perform comprehensive quality control:

• Filtering: Remove variants with call rate <95%, Hardy-Weinberg equilibrium p < 1×10⁻⁶, or excessive heterozygosity
• Annotation: Annotate variants using databases like gnomAD for population frequency, CADD for predicted deleteriousness, and REVEL for missense pathogenicity [36]
• Functional Prediction: Classify variants by functional impact (protein-truncating, missense, synonymous)

For the UK Biobank WGS analysis that identified IRS2 associations, researchers tested three variant categories: high-confidence protein-truncating variants, and missense variants with REVEL scores >0.5 or >0.7 [36].

Variant Filtering and Aggregation

Filter and aggregate rare variants within gene boundaries:

• Frequency Threshold: Apply minor allele frequency (MAF) filter of <0.1% for ultra-rare variants or <1% for rare variants [21]
• Functional Selection: Focus on putative high-impact variants (protein-truncating, deleterious missense) based on functional predictions
• Gene Aggregation: Collapse qualifying variants within each gene into a single burden score per individual

In the blood biomarker study, gene-based scores were calculated "based on the burden of rare functional variants and allele frequency" [17]. This approach increases power by aggregating signals from multiple rare variants within the same gene.

Association Testing

Perform statistical association between gene burden scores and phenotypes:

• Regression Models: Use linear regression for quantitative traits (e.g., biomarker levels) or logistic regression for binary traits (e.g., disease status)
• Covariate Adjustment: Include age, sex, genetic principal components, and relevant technical covariates
• Multiple Testing Correction: Apply Bonferroni correction based on the number of genes tested, with exome-wide significance threshold of p < 6.15×10⁻⁷ for 19,457 protein-coding genes [36]

For the UK Biobank analysis of BMI and T2D, the burden test framework examined 81,350 tests across different variant masks, requiring stringent multiple testing correction [36].

Validation and Sensitivity Analyses

• Replication: Test significant associations in independent cohorts (e.g., All of Us study, Icelandic populations)
• Sensitivity Analyses:
  • Leave-one-out: Verify no single variant drives the association
  • Common variant adjustment: Adjust for nearby common variants using polygenic risk scores
  • Ancestry-specific analyses: Evaluate consistency across ancestral backgrounds

The IRS2 discovery was replicated in the All of Us study, and sensitivity analyses confirmed the association was not driven by single variants or tagged by common variants [36].

Power Considerations for Study Design

The power of gene-based burden tests depends on several key factors:

Table 3: Factors Affecting Burden Test Power

Factor	Impact on Power	Recommendation
Sample Size	Critical determinant; >100,000 samples recommended for rare variants [21]	Leverage biobank-scale data (UK Biobank, All of Us)
Proportion of Causal Variants	Higher proportion increases power of burden tests [21]	Use functional annotations to enrich for causal variants
Variant Effect Sizes	Larger effects increase power; burden tests sensitive to effect direction	Focus on protein-truncating variants with large predicted effects
Variant Frequency Spectrum	Ultra-rare variants (MAF < 0.1%) often have larger effects	Use frequency thresholds appropriate to sample size

Analytical calculations and simulations based on 378,215 UK Biobank participants revealed that "aggregation tests are more powerful than single-variant tests only when a substantial proportion of variants are causal" [21]. For example, when protein-truncating variants, deleterious missense variants, and other missense variants have 80%, 50%, and 1% probabilities of being causal respectively, aggregation tests are more powerful for >55% of genes with n=100,000 and heritability of 0.1% [21].

The Scientist's Toolkit

Research Reagent Solutions

Table 4: Essential Research Reagents and Resources

Reagent/Resource	Function	Examples/Specifications
Whole Genome Sequencing	Comprehensive variant detection	Illumina NovaSeq, 30x coverage [3] [36]
Functional Annotation Databases	Variant prioritization	gnomAD (population frequency), CADD (deleteriousness), REVEL (missense pathogenicity) [36]
Burden Test Software	Statistical association testing	geneBurdenRD [3], COBT [5], SAIGE-GENE [17]
Biobank Datasets	Large-scale discovery cohorts	UK Biobank (∼500,000 participants), All of Us (∼1 million target) [36]
Phenotype Databases	Quantitative trait measurements	Blood biomarkers, electronic health records, medication data [17]
Stiripentol	Stiripentol, CAS:49763-96-4, MF:C14H18O3, MW:234.29 g/mol	Chemical Reagent
Sulforaphen	Sulforaphen, CAS:592-95-0, MF:C6H9NOS2, MW:175.3 g/mol	Chemical Reagent

Signaling Pathways Identified Through Burden Tests

Gene-based burden tests have illuminated several key pathways in cardiometabolic disease:

Cardiometabolic Signaling Pathways

Gene-based burden tests for rare variants have substantially advanced our understanding of cardiometabolic diseases and blood biomarker regulation. Through large-scale applications in biobank datasets, this approach has successfully identified novel disease genes, clarified biological mechanisms, and revealed potential therapeutic targets. The continued growth of genomic resources, combined with refined analytical methods and functional validation, promises to further enhance the impact of gene-based burden testing in cardiometabolic research and drug development.

The success stories presented here demonstrate that rare variants with large effects contribute significantly to the genetic architecture of complex cardiometabolic traits. As sample sizes increase and functional annotations improve, gene-based burden tests will remain an essential tool for uncovering the missing heritability of complex diseases and advancing precision medicine approaches in cardiometabolic health.

Serious adverse drug reactions (SADRs) represent a major clinical challenge, causing significant morbidity and mortality worldwide. It is estimated that SADRs occur in 6.2–6.7% of hospitalized patients in the United States, resulting in over 2 million incidents and more than 100,000 deaths annually [37]. These reactions are traditionally classified as Type A (pharmacologically predictable and dose-dependent) or Type B (idiosyncratic and unpredictable), with the latter proving particularly problematic for drug development and clinical management [37].

The integration of genetic evidence through burden testing of rare variants offers a transformative approach to predicting and preventing drug side effects. Gene-based burden tests address a fundamental limitation of traditional genome-wide association studies (GWAS), which primarily focus on common variants and have limited power for detecting associations with rare genetic variants [21]. By aggregating the effects of multiple rare variants within genes, burden tests can uncover genetic determinants of SADR susceptibility that would otherwise remain undetected [14].

The SE-GPS framework (Side Effect - Genomic Prediction Score) proposed in this protocol represents a novel methodology that leverages burden testing principles specifically for drug safety prediction. By integrating rare variant aggregation with pharmacogenomic insights, SE-GPS aims to identify patients at elevated risk for adverse drug events before treatment initiation, ultimately enabling more personalized and safer therapeutic interventions.

Technical Foundation: Principles of Gene-Based Burden Testing

Statistical Framework for Burden Tests

Gene-based burden tests operate on the principle that aggregating signals from multiple rare variants within a functional unit (typically a gene) can improve statistical power for detecting associations with traits or drug responses. The core methodology involves calculating a burden score for each individual by summing the minor allele counts across specified variants within a gene [21]. For an individual i with genotypes G~ij~ (coded as 0, 1, or 2 copies of the minor allele) at M variants in a gene, the burden score is computed as:

[ Burdeni = \sum{j=1}^{M} wj G{ij} ]

where w~j~ represents weights assigned to each variant based on functional impact or frequency [21]. The association between this burden score and a phenotype (such as drug response) is then tested using regression models.

The statistical power of burden tests depends critically on several factors: the proportion of causal variants among those aggregated, the effect sizes of these causal variants, the sample size, and the underlying genetic architecture [21]. Burden tests generally outperform single-variant tests when a substantial proportion of the aggregated variants are causal and exhibit effects in the same direction [21].

Comparison of Burden Tests with Single-Variant Approaches

The choice between burden tests and single-variant tests represents a key strategic decision in genetic association studies of drug side effects. Each approach has distinct advantages and limitations:

Table 1: Comparison of Single-Variant Tests vs. Burden Tests for SADR Identification

Feature	Single-Variant Tests	Burden Tests
Statistical Power	Higher for common variants with large effects	Higher for rare variants with convergent effects
Variant Scope	Tests each variant independently	Aggregates multiple variants within functional units
Genetic Architecture	Optimal for monogenic or oligogenic traits	Optimal for polygenic susceptibility with rare variants
Multiple Testing Burden	High (millions of tests)	Moderate (thousands of genes)
Interpretability	Direct variant-trait association	Gene-level association integrating multiple hits
Sample Size Requirements	Large samples needed for rare variants	More efficient for rare variants when appropriately aggregated

Single-variant tests maintain superiority for identifying associations with common variants or when only a small proportion of the aggregated rare variants are truly causal [21]. However, for SADRs where susceptibility may arise from the cumulative effect of multiple rare variants in genes involved in drug metabolism, transport, or target pathways, burden tests provide a more powerful framework [14].

SE-GPS Framework Implementation Protocol

Stage 1: Cohort Selection and Phenotyping

Objective: Assemble a well-characterized cohort with documented drug exposure and adverse event outcomes.

Procedures:

• Case Ascertainment: Identify individuals who experienced a specific SADR following drug exposure. Utilize standardized case definitions, such as:

  • Drug-Induced Liver Injury (DILI): Elevation of serum alanine aminotransferase (ALT) >5× upper limit of normal (ULN) or alkaline phosphatase >2× ULN, with appropriate temporal relationship to drug exposure [37].
  • Statin-Induced Myotoxicity: Creatine kinase elevation >10× ULN with muscle symptoms temporally related to statin therapy [37].
  • Drug-Induced Long QT Syndrome: QTc interval prolongation >500 ms or >60 ms increase from baseline following drug initiation [37].

• Control Selection: Identify controls with similar drug exposure but without adverse effects, matched for key covariates including age, sex, concomitant medications, and clinical comorbidities [37].

• Clinical Data Collection: Document detailed medication history (dose, duration, timing), clinical laboratory values, comorbidities, and demographic information using structured data collection forms.

• Sample Size Considerations: For rare variant analyses, larger sample sizes are critical. For SADRs with incidence of 1:1,000 to 1:10,000, sample sizes of several hundred cases and thousands of controls are typically required to achieve adequate power [14].

Stage 2: Genomic Data Processing and Quality Control

Objective: Generate high-quality genetic data suitable for burden testing analyses.

Procedures:

• Sequencing Approach: Perform whole-exome or whole-genome sequencing using established platforms (e.g., Illumina NovaSeq). Maintain consistent sequencing depth (>30× mean coverage) across all samples.

• Variant Calling: Implement standardized variant calling pipelines such as the GATK Best Practices workflow [14]:

  • Alignment to reference genome (GRCh38) using BWA-MEM
  • Base quality score recalibration and indel realignment
  • Variant calling with GATK HaplotypeCaller
  • Joint genotyping across all samples

• Quality Control Filters:

  • Sample-level: Exclude samples with call rate <95%, contamination >3%, or sex discrepancies
  • Variant-level: Retain variants with call rate >95%, QUAL >30, and PASS filter status
  • Relatedness: Identify and remove related individuals (IBD >0.1875)

• Variant Annotation: Annotate variants using Variant Effect Predictor (VEP) with additional pathogenicity predictions from PolyPhen-2, SIFT, and CADD [14].

Stage 3: Variant Selection and Mask Definition

Objective: Define biologically informed variant masks for burden testing.

Procedures:

• Variant Category Definition: Classify variants based on functional impact and frequency:

  • Protein-truncating variants (PTVs): Stop-gain, frameshift, splice-site variants
  • Deleterious missense: Missense variants with CADD score >20 or REVEL score >0.5
  • Other missense: Remaining missense variants not meeting deleterious criteria

• Frequency Thresholds: Apply population-specific frequency filters based on gnomAD or other reference databases:

  • Rare variants: MAF <0.1% in any population
  • Ultra-rare variants: MAF <0.01% in any population

• Mask Specification: Create predefined variant masks for testing:

  • PTV-only mask
  • PTV + deleterious missense mask
  • Gene-specific masks informed by biological knowledge (e.g., metabolizing enzyme genes for specific drugs)

• Variant Weighting Schemes: Implement frequency-dependent weighting such as:

  • Beta(MAF, 1, 25) weights as used in SKAT
  • Madsen-Browning weights
  • Functional prediction-based weights (e.g., CADD-scaled)

Stage 4: Burden Test Implementation and Analysis

Objective: Perform gene-based burden tests to identify genetic associations with SADRs.

Procedures:

• Test Selection: Implement multiple burden tests to accommodate different genetic architectures:

  • Collapsing burden test: Equal weights for all qualifying variants
  • Weighted burden test: Frequency-weighted or function-weighted variants
  • SKAT: Variance-component test robust to mixed effect directions
  • SKAT-O: Optimally combines burden and SKAT tests

• Software Implementation: Utilize established software packages:

  • TRAPD (Test Rare vAriants with Public Data) for case-control burden testing [14]
  • OWC (Optimally Weighted Combination) for leveraging multiple weighting schemes [38]
  • STAAR for comprehensive variant set tests

• Control Population: Use public databases (gnomAD, UK Biobank) as controls when internal controls are limited, implementing stringent quality control to address batch effects and platform differences [14].

• Covariate Adjustment: Adjust for potential confounding factors including genetic ancestry (principal components), age, sex, and relevant clinical covariates in regression models.

• Significance Thresholding: Apply gene-based multiple testing correction (Bonferroni correction for ~20,000 genes: p<2.5×10^-6^) or false discovery rate control (FDR <0.05).

Stage 5: SE-GPS Risk Score Development

Objective: Translate burden test results into clinically actionable risk scores.

Procedures:

• Multi-gene Integration: Combine signals from multiple associated genes into a polygenic risk score: [ SE-GPSi = \sum{g=1}^{G} \betag \times Burden{ig} ] where (\betag) is the effect size estimate for gene *g* and (Burden{ig}) is the burden score for individual i in gene g.

• Model Validation: Perform cross-validation or bootstrap validation to assess model performance and prevent overfitting.

• Performance Metrics: Evaluate discrimination (AUC-ROC, C-statistic), calibration (Hosmer-Lemeshow test), and reclassification (NRI) metrics.

• Clinical Implementation: Establish risk strata based on score distribution (e.g., low, intermediate, high risk) with corresponding clinical management recommendations.

Application to Specific Adverse Drug Reactions

Drug-Induced Liver Injury (DILI)

Drug-induced liver injury represents a major challenge in drug development and clinical practice, with genetic factors playing a significant role in susceptibility. The SE-GPS framework can be applied to DILI through the following specific approaches:

Candidate Genes: Focus on genes involved in drug metabolism (CYPs, UGTs), transport (ABCB1, ABCC2), and immune response (HLA genes) [37]. For example, polymorphisms in UGT1A1 have been associated with tranilast-induced hyperbilirubinemia, and CYP2E1 and NAT2 variants with isoniazid-induced hepatotoxicity [37].

Variant Mask: Create a DILI-specific mask incorporating:

• Protein-truncating variants in drug metabolism genes
• Deleterious missense variants in HLA genes (particularly DRB1*07:01 associated with ximelagatran hepatotoxicity) [37]
• Rare variants in bile salt export pump (BSEP/ABCB11) and other hepatobiliary transporters

Validation: Replicate associations in independent DILI cohorts (e.g., DILI Network, International DILI Consortium) and perform functional validation in hepatocyte models.

Statin-Induced Myotoxicity

Statin-associated muscle symptoms represent one of the most common reasons for discontinuation of therapy, with genetic factors influencing susceptibility.

Candidate Genes: Prioritize genes involved in statin pharmacokinetics (SLCO1B1, ABCG2) and pharmacodynamics (HMGCR, COQ2), as well as muscle-specific pathways.

Variant Mask: Develop a myotoxicity-specific mask including:

• Rare deleterious variants in SLCO1B1 beyond the common c.521T>C polymorphism
• Protein-truncating variants in genes involved in coenzyme Q10 biosynthesis
• Rare variants in skeletal muscle ion channels and mitochondrial function genes

Clinical Application: Combine the SE-GPS score with clinical risk factors (age, diabetes, renal function) for comprehensive risk assessment.

Table 2: Genetic Risk Factors for Selected Serious Adverse Drug Reactions

Adverse Drug Reaction	Implicated Genes	Variant Types	Biological Mechanism	Clinical Implications
Drug-Induced Liver Injury	HLA-DRB1*07:01, UGT1A1, CYP2E1, NAT2	Missense, regulatory	Altered drug metabolism, immune activation	Avoid ximelagatran in HLA-DRB1*07:01 carriers; dose adjustments for isoniazid in slow acetylators
Statin-Induced Myotoxicity	SLCO1B1, ABCG2, COQ2, HMGCR	Nonsynonymous, loss-of-function	Reduced hepatic uptake, mitochondrial dysfunction	Alternative statins or reduced doses in SLCO1B1 risk variant carriers
Drug-Induced Long QT Syndrome	KCNH2, KCNQ1, KCNE1, SCNSA	Missense, splice-site	Impaired cardiac repotassium current	Avoid QT-prolonging drugs in variant carriers; enhanced monitoring
Abacavir Hypersensitivity	HLA-B*57:01	Coding, regulatory	Immune recognition of drug-peptide complex	Pretreatment screening eliminates hypersensitivity risk

Cross-Ancestry Considerations in SE-GPS Implementation

Genetic ancestry significantly influences pharmacogenomic risk profiles, with distinct geographical patterns observed for drug toxicity risk [39]. Admixed American and European populations demonstrate higher overall risk proximity for drug toxicity, while East Asian and Oceanian populations display relatively protective genetic profiles [39].

Ancestry-Specific Adjustments:

• Develop ancestry-specific variant frequency thresholds based on corresponding gnomAD populations
• Account for differences in LD structure and haplotype background across populations
• Consider population-specific variant effects when supported by sufficient data
• Validate SE-GPS models in multiple ancestry groups before broad clinical implementation

Table 3: Key Research Reagents and Computational Tools for SE-GPS Implementation

Resource Category	Specific Tools/Databases	Primary Function	Key Features
Variant Calling	GATK, BCFtools, BWA	Sequence alignment, variant discovery	Base quality recalibration, indel realignment, joint genotyping [14]
Variant Annotation	VEP, ANNOVAR, SnpEff	Functional consequence prediction	Gene context, protein effect, regulatory regions [14]
Pathogenicity Prediction	PolyPhen-2, SIFT, CADD	Deleteriousness assessment	Evolutionary conservation, structural impact [14]
Variant Frequency Databases	gnomAD, TOPMed, UK Biobank	Population frequency data	Ancestry-specific allele frequencies, quality metrics [14]
Burden Testing Software	TRAPD, OWC, STAAR, SKAT/SKAT-O	Gene-based association tests	Public control integration, multiple weighting schemes [14] [38]
Clinical Interpretation	PharmGKB, ClinVar, DrugBank	Pharmacogenomic knowledge	Drug-gene interactions, clinical guidelines [37]
Visualization	IGV, LocusZoom, R/bioconductor	Data exploration and presentation	Genomic context, association signals, quality metrics

Integration Pathways from Genetic Discovery to Clinical Application

The translation of burden test findings into clinically actionable SE-GPS implementations requires a structured pathway from genetic discovery to clinical integration:

• Genetic Discovery Phase: Initial burden tests identify genes associated with SADR susceptibility across multiple variant masks and statistical approaches.

• Functional Validation: Confirm biological mechanisms through:

  • In vitro assays of protein function for identified rare variants
  • Drug transport/metabolism studies in cellular models
  • Animal models replicating key genetic perturbations

• Risk Model Development: Refine SE-GPS algorithms through:

  • Integration of multiple genetic signals into composite scores
  • Determination of optimal risk thresholds for clinical decision-making
  • Development of clinical reporting formats

• Clinical Validation: Assess real-world performance in:

  • Prospective clinical trials evaluating SE-GPS-guided therapy
  • Healthcare systems implementing pre-prescription genetic screening
  • Diverse patient populations to ensure generalizability

• Guideline Development and Implementation: Incorporate SE-GPS into:

  • Drug labeling and package inserts
  • Professional society guidelines for high-risk medications
  • Electronic health record decision support tools

This structured approach ensures that genetic discoveries arising from burden tests are rigorously evaluated and effectively implemented to improve medication safety.

The SE-GPS framework represents a powerful methodology for predicting adverse drug reactions by integrating gene-based burden tests of rare variants with pharmacogenomic principles. As genetic datasets continue to expand in size and diversity, and as functional annotation of rare variants improves, the accuracy and clinical utility of SE-GPS predictions will continue to advance.

Future developments in this field will likely include the integration of transcriptomic and proteomic data to refine risk predictions, the application of machine learning methods to identify complex interaction effects, and the expansion of SE-GPS to diverse ancestral populations to ensure equitable benefits across all patient groups.

By providing a standardized protocol for implementing burden tests in drug safety assessment, this SE-GPS framework enables researchers and drug development professionals to systematically identify genetic risk factors for adverse drug reactions and translate these findings into clinically actionable strategies for personalized medicine.

Navigating Analytical Pitfalls: Ensuring Robust and Reproducible Results

Addressing Population Stratification and Batch Effects

In gene-based burden tests for rare variants, controlling for technical and biological confounders is paramount for ensuring the validity of association findings. Population stratification—systematic ancestry differences between case and control groups—and batch effects—non-biological variations introduced by experimental procedures—can both induce spurious associations if not properly addressed [40]. These issues are particularly acute in rare variant analysis due to the low frequency of the genetic variants of interest, which can be correlated with ancestry or batch membership. This Application Note details the protocols and methodologies researchers can employ to mitigate these threats, framed within the context of a broader thesis on rare variant research.

Quantitative Comparison of Correction Methods

The choice of statistical method to correct for population stratification is critical and depends on factors such as sample size and the complexity of the population structure. The table below summarizes the type I error rates (the false positive rate) of different correction methods under various stratification scenarios, as evaluated through a large simulation study using real exome data [40].

Table 1: Empirical Type I Error Rates (x10⁻⁶) at α=2.5x10⁻⁶ for Various Correction Methods under Different Stratification Scenarios [40]

Sample Size (Cases:Controls)	Stratification Scenario	Uncorrected Test	Principal Components (PC)	Linear Mixed Models (LMM)	Local Permutation (LocPerm)
750:4250 (15% cases)	Within-Continent (Moderate)	180.0	3.8	2.5	2.5
750:4250 (15% cases)	Between-Continent (Severe)	450.0	350.0	4.5	3.2
50:283 (15% cases)	Within-Continent (Moderate)	175.0	15.8 (with 100 controls)	12.5 (with 1000 controls)	3.5
50:283 (15% cases)	Between-Continent (Severe)	480.0	380.0 (with 100 controls)	10.5 (with 1000 controls)	3.8

Key Interpretation: The data shows that without correction, type I error is massively inflated. While PC and LMM methods can be effective, their performance varies with sample size and stratification severity. The Local Permutation (LocPerm) method maintains robust error control across all scenarios, including those with small case numbers, which are common in studies of rare diseases [40].

Detailed Experimental Protocols

Protocol 1: Controlling Population Stratification in Cohort-Level Analysis

This protocol is designed for a single study cohort where individual-level genetic data is available [40].

• Quality Control (QC) and Genotype Processing:
  • Perform standard QC on whole-exome or whole-genome sequencing data. Retain only coding variants with a depth of coverage (DP) > 8, genotype quality (GQ) > 20, minor read ratio (MRR) > 0.2, and a call rate > 95%.
  • Exclude related individuals based on kinship coefficients (e.g., King's kinship > 0.1875) to ensure a dataset of unrelated samples.
• Population Structure Assessment:
  • Merge the study cohort with reference populations (e.g., 1000 Genomes) to provide a global ancestry context.
  • Perform a Principal Component Analysis (PCA) on common variants (MAF > 1%) to capture the major axes of genetic variation.
  • Visually inspect the PCA plot (PC1 vs. PC2) to identify population outliers and major ancestry clusters.
• Association Testing with Covariate Adjustment:
  • For burden tests (e.g., CAST test), include the top principal components (typically 5-10) as covariates in the regression model to adjust for broad-scale population stratification.
  • Alternative Method: Apply a Linear Mixed Model (LMM), which uses a genetic relatedness matrix (GRM) to account for finer-scale relatedness and population structure simultaneously [18].
• Validation with Local Permutation (LocPerm):
  • For studies with small sample sizes or complex stratification where PC adjustment may be insufficient, employ the LocPerm method.
  • This method involves permuting case-control labels within genetically similar, locally defined neighborhoods, thereby preserving the local population structure while breaking the association with the phenotype. This provides a robust null distribution for significance testing [40].

This protocol is for combining results from multiple cohorts using summary statistics, a common practice in consortium studies. The Meta-SAIGE and REMETA methods are highlighted for their efficiency and accurate error control [18] [41].

• Preparing Summary Statistics and LD Matrices (Per Cohort):
  • Single-Variant Association: Use a logistic or linear mixed model (e.g., SAIGE, REGENIE) that accounts for sample relatedness to generate per-variant score statistics (S) and their variances for each trait [18] [41].
  • LD Matrix Calculation: Compute a sparse linkage disequilibrium (LD) matrix (Ω) for each gene or region. This matrix contains the pairwise correlations between variants.
  • Key Innovation: Unlike other methods, Meta-SAIGE and REMETA use a single, sparse LD matrix per study that is not phenotype-specific. This matrix can be re-used across all phenotypes, drastically reducing computational cost and storage requirements [18] [41]. The formula for storage requirements demonstrates this efficiency: O(MFK + MKP) for Meta-SAIGE/REMETA vs. O(MFKP + MKP) for methods requiring trait-specific LD matrices, where M is variants, F is non-zero LD pairs, K is cohorts, and P is phenotypes [18].
• Combining Summary Statistics (Meta-Analysis):
  • Consolidate score statistics from all cohorts. For binary traits with case-control imbalance, apply statistical corrections like the Saddlepoint Approximation (SPA) and genotype-count-based SPA (GC-SPA) to accurately estimate the null distribution and control type I error [18].
  • The covariance matrix of the combined score statistics is calculated as Cov(S) = V¹ᐟ² Cor(G) V¹ᐟ², where Cor(G) is the correlation matrix from the sparse LD matrix (Ω), and V is the diagonal matrix of the variance of S [18].
• Gene-Based Burden Testing:
  • With the combined per-variant statistics and covariance matrix, perform gene-based burden tests, SKAT, or SKAT-O.
  • Handle ultra-rare variants (e.g., MAC < 10) by collapsing them to enhance power and error control.
  • Combine P values from different functional annotations and MAF cutoffs using an omnibus test like the Cauchy combination method [18].

Protocol 3: Correcting for Batch Effects in Microbiome Studies (Analogous Workflow)

While from a different field, the principles of batch effect correction in microbiome data are analogous to genetic data. The following protocol uses a composite quantile regression approach [42].

• Model Systematic Batch Effects:
  • For operational taxonomic unit (OTU) count data Yijg for sample i, OTU j, and batch g, fit a negative binomial regression model: log(μijg) = σj + Xiβj + γjg + logNi where μijg is the expected count, σj is the OTU-specific baseline, Xi are biological covariates, βj are their coefficients, γjg is the mean batch effect, and Ni is the library size [42].
  • Adjust the counts by subtracting the estimated batch effect: log(μ*ij) = log(μijg) - γjg [42].
• Address Non-Systematic Batch Effects:
  • Apply composite quantile regression to adjust the distribution of OTUs in each batch to match a reference batch. This step handles batch effects that vary in strength depending on the specific OTU and sample.
  • Select a reference batch using a method like the Kruskal-Wallis test to identify the most representative batch.
• Validation:
  • Evaluate the success of batch correction using Principal Coordinates Analysis (PCoA) plots and metrics like the Average Silhouette Coefficient, which should show batches mixing well after correction [42].

Workflow Visualization

The diagram below illustrates the integrated workflow for meta-analysis of gene-based burden tests while controlling for population stratification and batch effects, as described in the protocols.

The Scientist's Toolkit: Research Reagent Solutions

The following table lists key software and methodological tools essential for implementing the protocols described above.

Table 2: Essential Tools for Addressing Confounders in Rare Variant Studies

Tool Name	Type	Primary Function	Key Feature
SAIGE/SAIGE-GENE+ [18]	Software	Association testing for binary traits	Uses Saddlepoint Approximation (SPA) to control for case-control imbalance and relatedness.
Meta-SAIGE [18]	Software	Rare variant meta-analysis	Extends SAIGE-GENE+; uses a single LD matrix per cohort for all phenotypes.
REMETA [41]	Software	Gene-based meta-analysis	Uses rescaled reference LD matrices; integrates with REGENIE for efficient testing.
LocPerm [40]	Method	Type I error control	Local permutation for robust error control in small samples or complex stratification.
Composite Quantile Regression [42]	Method	Batch effect correction	Corrects both systematic and non-systematic batch effects in count data (e.g., microbiome).
Reference LD Matrix [18] [41]	Data Resource	Covariance for summary statistics	Pre-computed, sparse LD matrix for a cohort, reusable across phenotypes to save storage.
Saddlepoint Approximation (SPA) [18]	Statistical Method	P-value calculation	Provides accurate P-values for score tests under severe case-control imbalance.
Tanshinone I	Tanshinone I		Bench Chemicals
Syringin	Syringin (Eleutheroside B) - CAS 118-34-3 - For Research	High-purity Syringin for research. Explore its applications in cardiovascular, metabolic, and neuroscience studies. For Research Use Only. Not for human use.	Bench Chemicals

In the field of genetic association studies, gene-based burden tests for rare variants have become a cornerstone for identifying associations with complex diseases. However, a significant statistical challenge arises when these tests are applied to binary traits with extremely unbalanced case-control ratios, a common scenario in studies of rare diseases or low-prevalence conditions within large biobanks. Traditional association tests, such as the score test that relies on normal approximation, often fail under these conditions, leading to severely inflated Type I error rates [43] [44]. This inflation occurs because the underlying distribution of the test statistic is highly skewed in unbalanced designs, and the normal approximation, which depends only on the first two moments (mean and variance), cannot accurately capture this skewness, especially in the tails of the distribution where p-values are calculated [43].

The Saddlepoint Approximation (SPA) has emerged as a superior method to overcome this limitation. Contrary to normal approximation, SPA uses the entire cumulant-generating function (CGF), effectively incorporating higher moments such as skewness into its estimation of the null distribution [43] [44]. This methodological advance allows SPA to accurately calibrate p-values even when case-control ratios are as extreme as 1:600 [43]. The practical implication is profound: researchers can now conduct large-scale phenome-wide association studies (PheWAS) across thousands of binary traits in biobanks with confidence, knowing that the false positive rate is well-controlled. This reliability is crucial for both the discovery of novel gene-trait associations and the subsequent translation of these findings into drug development pipelines.

Methodological Foundation: From Normal Approximation to SPA

The Pitfalls of Standard Methods

The inadequacy of standard methods in unbalanced designs is not trivial. The logistic mixed model, implemented in tests like the Generalized Mixed Model Association Test (GMMAT), assumes that score test statistics asymptotically follow a Gaussian distribution. However, this assumption is invalidated by unbalanced case-control ratios, leading to poor calibration [44]. Similarly, linear mixed models (LMMs), while widely used for their computational efficiency, are not designed for binary traits and exhibit large Type I error rates in such scenarios [44] [45]. This creates a pressing need for a method that is both computationally feasible for biobank-scale data and statistically robust for unbalanced binary traits.

The Saddlepoint Approximation (SPA) Explained

SPA addresses the core of the problem by providing a much more accurate approximation of the null distribution. The key advantage is its relative error bound of O(n^{-3/2}), a substantial improvement over the O(n^{-1/2}) convergence rate of normal approximation [43]. In practice, this means that SPA provides highly accurate p-value calculations far into the tails of the distribution, which is exactly where significance testing occurs.

The implementation typically involves calculating a score statistic (S) for each genetic variant. The core of SPA lies in using the estimated CGF of this statistic, K(t), and its derivatives to approximate the probability of observing a score statistic as extreme as the one calculated [43]. The distribution of the score statistic S at a value s can be approximated as: Pr(S < s) ≈ F~(s) = Φ{w + (1/w) * log(v/w)} where w and v are derived from the saddlepoint equation and Φ is the standard normal cumulative distribution function [43]. This formula, while complex, is the engine that enables accurate error control.

Table 1: Comparison of Approximation Methods for Association Testing

Feature	Normal Approximation	Saddlepoint Approximation (SPA)
Statistical Basis	Uses only mean and variance (first two moments)	Uses the entire cumulant-generating function (all moments)
Error Bound	O(n^{-1/2})	O(n^{-3/2})
Performance in Unbalanced Designs	Poor; highly inflated Type I error	Excellent; well-controlled Type I error
Computational Cost	Low	Higher than normal approximation, but optimized versions (e.g., fastSPA) are efficient
Suitable for	Balanced studies, common variants	Unbalanced case-control studies, rare and low-frequency variants

Protocol: Implementing SPA in Genetic Association Studies

The SAIGE Workflow: A Step-by-Step Guide

The Scalable and Accurate Implementation of GEneralized mixed model (SAIGE) is a robust software tool that effectively implements SPA for large-scale genetic association studies [44] [45]. Its protocol is divided into two main steps, designed for efficiency and accuracy with large sample sizes.

Step 1: Fitting the Null Logistic Mixed Model

• Objective: To estimate the variance components and other parameters of the null model (without genetic variants) that account for sample relatedness and population structure.
• Procedure:
  • Input: Genotype matrix, phenotypic data (binary trait), and covariates (e.g., age, sex, principal components).
  • Model Fitting: A null logistic mixed model is fit using the Average Information Restricted Maximum Likelihood (AI-REML) algorithm.
  • Computational Optimization: To handle hundreds of thousands of samples, SAIGE employs the preconditioned conjugate gradient (PCG) method to solve linear systems iteratively without needing to store the massive Genetic Relationship Matrix (GRM) in memory [44]. The raw genotypes are stored in a compact binary format, drastically reducing memory usage (e.g., from ~669 GB to ~9.5 GB for N=400,000) [44].
  • Output: Estimated variance components, model parameters, and predicted random effects for each individual.

Step 2: Association Testing with SPA Calibration

• Objective: To test each genetic variant for association with the binary trait, using SPA to calibrate the p-values accurately.
• Procedure:
  • Variance Ratio Calibration: A constant variance ratio is estimated from a subset of variants and used to adjust the variance of the score statistics for all variants, which helps account for the polygenic background [44].
  • Score Test & SPA Application: For each variant, a score statistic is computed. SAIGE then applies SPA (specifically, a faster version similar to fastSPA [43]) to approximate the null distribution and calculate an accurate p-value [44]. The fastSPA algorithm is particularly efficient for rare variants as its computation cost depends on the number of minor allele carriers rather than the total sample size [43].
  • Output: Association p-values and effect size estimates for all tested variants.

The following workflow diagram illustrates the key stages of the SAIGE protocol:

Application to Rare Variant Meta-Analysis: Meta-SAIGE

The principles of SPA have been successfully extended to the meta-analysis of gene-based rare variant tests through the Meta-SAIGE method [18]. This is critical because meta-analysis enhances power by combining data from multiple cohorts, but existing methods like MetaSTAAR can exhibit inflated Type I error for low-prevalence traits [18].

The Meta-SAIGE protocol involves:

• Preparation: Each cohort uses SAIGE to generate per-variant score statistics and a sparse linkage disequilibrium (LD) matrix. The LD matrix is not phenotype-specific and can be reused across different phenotypes, drastically reducing computational and storage burdens [18].
• Summary Statistics Combination: Score statistics from all cohorts are combined. To further improve error control, Meta-SAIGE applies a genotype-count-based SPA in addition to the SPA applied within each cohort [18].
• Gene-Based Testing: Burden, SKAT, and SKAT-O tests are performed on the combined summary statistics, similar to the SAIGE-GENE+ method for individual-level data.

Table 2: Benchmarking SAIGE Performance in the UK Biobank

Metric	Value for Coronary Artery Disease (N=408,458)	Context and Comparison
Projected Computation Time	517 CPU hours	For 71 million variants; faster than GMMAT, slightly slower than BOLT-LMM [44]
Memory Usage	10.3 Gb	Feasible for large-scale analysis; lower than GMMAT and GEMMA [44]
Type I Error Control	Accurate	Controls inflation for traits with extreme case-control imbalance (e.g., <1:100) [44] [18]
Power	High	Comparable to joint analysis of individual-level data when used in meta-analysis (Meta-SAIGE) [18]

To implement the protocols described, researchers require access to specific computational tools and reference data. The following table details key components of the research toolkit.

Table 3: Research Reagent Solutions for SPA-Based Association Studies

Tool/Resource	Type	Primary Function	Key Feature
SAIGE [44] [45]	Software Tool	Scalable association testing for binary traits	Integrates logistic mixed models with SPA calibration for unbalanced case-control studies
SPAtest [43]	R Package	Score test with SPA for binary phenotypes	Provides the fastSPA algorithm for efficient computation with rare variants
Meta-SAIGE [18]	Software Tool	Rare variant meta-analysis	Extends SPA to meta-analysis, controlling Type I error across cohorts
REGENIE/REMETA [41]	Software Tool	Whole-genome regression & meta-analysis	Uses a reference LD matrix for efficient gene-based tests from summary statistics
Population Reference Data (e.g., gnomAD) [46]	Data Resource	Aggregate variant frequencies	Provides population-level allele frequencies for calibrating case-only tests and interpreting results

The integration of Saddlepoint Approximation into the statistical toolkit for genetic association studies represents a critical advancement. By moving beyond the limitations of normal approximation, SPA-based methods like SAIGE and Meta-SAIGE provide the robustness required to explore the genetic underpinnings of thousands of binary phenotypes in biobanks confidently. This robustness, combined with computational optimizations that make analysis of hundreds of thousands of samples feasible, empowers researchers and drug development professionals to discover and validate genetic associations with greater accuracy, ultimately accelerating the translation of genetic insights into therapeutic strategies.

Optimizing Power through Variant Quality Filtering and Functional Annotation

Gene-based burden tests have become a cornerstone in rare variant (RV) association studies for complex diseases and Mendelian disorders. These tests aggregate rare variants within a gene or functional unit to increase statistical power for detecting associations that would be elusive at the single-variant level [47] [14]. However, the effectiveness of burden testing is highly dependent on two critical factors: the quality of variant calling and the strategic incorporation of functional annotations. Variant quality filtering ensures that true biological signals are distinguished from technical artifacts, while functional annotation prioritizes variants with potentially deleterious biological consequences. This protocol details a comprehensive framework that integrates rigorous quality control measures with functionally-informed analytical methods to optimize the power of rare variant association studies, with particular emphasis on the gruyere (Genome-wide Rare Variant EnRichment Evaluation) method [47] [48].

Variant Quality Control Framework

Genotype-Level Quality Filtering

Sequencing artifacts and genotyping errors significantly reduce power in rare variant association tests, particularly for very low-frequency variants. Stringent genotype-level filtering is essential before proceeding to association testing. The following thresholds, empirically derived from whole-exome sequencing (WES) studies, substantially improve genotype concordance and transition/transversion (Ti/Tv) ratios [49].

Table 1: Recommended Genotype-Level Quality Filters

Filtering Metric	Threshold	Rationale
Depth (DP)	≥ 10 reads	Ensures sufficient coverage for accurate genotype calling
Genotype Quality (GQ)	≥ 20	Phred-scaled confidence > 99% in genotype call
Missingness (Call Rate)	> 95% per variant	Removes variants with excessive missing data

Implementation of these filters in a large Norwegian cohort (n=920) improved non-reference genotype concordance from 99.33% to 99.77% and the Ti/Tv ratio from 2.63 to 2.75, indicating substantial reduction of false positives [49]. These filters are particularly crucial for rare variants, where heterozygote-to-homozygote errors dramatically reduce statistical power in burden tests [49].

Variant-Level Quality Filtering

Following genotype filtering, variant-level filters address systematic biases and batch effects. The GATK Variant Quality Score Recalibration (VQSR) tool provides a foundation, but supplemental filtering is recommended for optimal results [49].

Table 2: Supplemental Variant-Level Quality Filters

Filtering Metric	Threshold	Application Note
Hardy-Weinberg Equilibrium (HWE)	p > 1 × 10-6	Identifies gross genotyping errors and population stratification
Batch Effects	Separate processing by capture kit/chemistry	Managed by analyzing batches separately; recovered 40.9% more high-quality variants in one study [49]

When using public databases (e.g., gnomAD) as controls, additional calibration is essential. The TRAPD (Test Rare vAriants with Public Data) method recommends using synonymous variants as internal controls to calibrate for ancestry differences, sequencing platform variations, and bioinformatic processing discrepancies [14].

Functional Annotation Strategies

Functional annotations provide biological context for prioritizing rare variants. Integrating multiple annotation sources significantly improves the identification of pathogenic variants.

Table 3: Key Functional Annotation Resources for Rare Variants

Annotation Category	Tool/Resource	Application in Burden Testing
Variant Effect Prediction	Variant Effect Predictor (VEP) [50]	Standardized annotation of consequence types (e.g., missense, LoF)
Deleteriousness Meta-Scores	CADD, Eigen [51]	Integrates multiple genomic features into a single measure of variant deleteriousness
Splicing Disruption	SpliceAI [47] [48]	Deep learning-based prediction of splicing alterations
Non-Coding Functional Effects	Enformer [47] [48]	Predicts cell-type-specific effects on transcription factor binding and chromatin state
Enhancer-Gene Links	Activity-by-Contact (ABC) Model [47] [48]	Maps non-coding variants to target genes using chromatin state and conformation data

The Gruyere Framework for Annotation Integration

The gruyere method provides a Bayesian probabilistic framework that moves beyond simple variant filtering to learn trait-specific weights for functional annotations genome-wide [47] [48]. This approach addresses a key limitation of standard burden tests, which typically apply static, one-size-fits-all annotation filters.

In the gruyere model, the genetic risk for an individual ( i ) from variants in gene ( g ) is modeled as:

[ \text{logit}(\mu{ig}) = Xi\alphag + G{ig}\beta_{gj} ]

where ( \mu{ig} ) represents the probability of disease, ( Xi ) are covariates, ( G{ig} ) are genotypes, ( \alphag ) are covariate weights, and ( \beta_{gj} ) are variant effects [47].

The key innovation is that variant effects are parameterized as:

[ \beta{gj} = wg \times wj \times \left( \tau0 + \sum{k=1}^{q} Z{gjk}\tau_k \right) ]

where:

• ( w_g ) = gene importance weight (learned)
• ( w_j ) = variant weight based on minor allele frequency
• ( Z_{gjk} ) = functional annotations for each variant
• ( \tau_k ) = annotation importance weights (learned globally across genes) [47]

This model enables gruyere to simultaneously discover associated genes while learning which functional annotations are most relevant for the specific trait being studied.

Integrated Protocols

Protocol 1: Quality Control for Burden Testing

Objective: Generate high-quality variant calls optimized for rare variant burden tests.

• Variant Calling: Process sequencing data through standard GATK Best Practices, including base quality score recalibration, indel realignment, and variant quality score recalibration (VQSR) [14] [49].
• Apply Genotype Filters: Filter genotypes using thresholds in Table 1 (DP ≥ 10, GQ ≥ 20). This can be implemented using bcftools filter or custom scripts.
• Apply Variant Filters: Remove variants failing HWE (p < 1 × 10^-6) and with call rate < 95%. Analyze samples from different capture kits or sequencing batches separately to manage batch effects [49].
• Control Cohort Calibration (if using public databases): Use the TRAPD framework [14]:
  • Annotate variants in cases and public controls (e.g., gnomAD) consistently.
  • Use synonymous variants to assess and calibrate technical differences between case and control datasets.
  • Apply adaptable variant quality filters based on this calibration.

Protocol 2: Functionally-Informed Burden Testing with Gruyere

Objective: Perform genome-wide burden testing integrating functional annotations to improve power for gene discovery.

• Define Testing Units:

  • For coding variants: Use canonical gene boundaries from GENCODE or RefSeq.
  • For non-coding variants: Use the Activity-by-Contact (ABC) model to define enhancer-gene links in relevant cell types (e.g., microglia for Alzheimer's disease) [47] [48].

• Annotate Variants: Generate a comprehensive annotation matrix (( Z )) for all rare variants (MAF < 0.01) including:

  • Predicted deleteriousness scores (CADD)
  • Splicing disruption (SpliceAI)
  • Cell-type-specific regulatory effects (Enformer)
  • Chromatin state predictions [47] [48]

• Initial Gene Selection: Perform an initial burden test (e.g., Functional Score Test) with a lenient significance threshold (nominal p < 0.05) to select candidate genes for gruyere analysis [47].

• Implement Gruyere Model:

  • Split data into training (80%) and test (20%) sets
  • Using the training set, fit the hierarchical Bayesian model to estimate:
    • Gene importance weights (( wg ))
    • Global functional annotation weights (( \tauk ))
    • Covariate coefficients (( \alpha_g ))
  • Evaluate model performance on the test set to ensure robust generalization [47]

• Interpret Results: Prioritize genes with significant ( wg ) values and examine the ( \tauk ) weights to identify functional annotations most predictive of disease risk.

The Scientist's Toolkit

Table 4: Essential Research Reagents and Computational Tools

Resource	Type	Function	Access
GATK [49]	Software Suite	Industry standard for variant discovery in high-throughput sequencing data	https://gatk.broadinstitute.org
gnomAD [14]	Control Database	Publicly available population frequency reference for variant filtering and calibration	https://gnomad.broadinstitute.org
VEP [50]	Annotation Tool	Predicts functional consequences of variants on genes, transcripts, and protein sequence	https://www.ensembl.org/info/docs/tools/vep/index.html
SpliceAI [47] [48]	Deep Learning Model	Quantifies variant effects on mRNA splicing using a deep residual neural network	https://github.com/Illumina/SpliceAI
Enformer [47] [48]	Deep Learning Model	Predicts variant effects on cell-type-specific gene regulation from DNA sequence	https://github.com/google-deepmind/deepmind-research/tree/master/enformer
ABC Model [47] [48]	Algorithm	Predicts enhancer-gene connections using chromatin state and conformation data	https://github.com/broadinstitute/ABC-Enhancer-Gene-Prediction
gruyere [47] [48]	R Package	Bayesian framework for rare variant association testing with functional annotations	https://github.com/daklab/gruyere
TRAPD [14]	Software Package	Implements burden testing of rare variants using public control databases	https://github.com/ThomasYeoLab/TRAPD

Application Notes

The integrated quality filtering and functional annotation approach has demonstrated significant utility across multiple domains:

• Alzheimer's Disease: Application of gruyere to whole-genome sequencing data from the Alzheimer's Disease Sequencing Project (7,966 cases, 13,412 controls) identified 15 significant genetic associations not detected by other RV methods, with deep learning-based annotations for splicing, transcription factor binding, and chromatin state being highly predictive of functional non-coding rare variants [47] [48].

• Cancer Genetics: Gene-based burden tests aggregating rare missense and loss-of-function variants across 130,991 cancer cases and 733,486 controls identified six novel cancer susceptibility genes, including BIK (prostate cancer), ATG12 (colorectal cancer), and PPP1R15A (breast cancer) [50] [10].

• Blood Biomarkers: Integration of gene-based burden scores with polygenic risk scores for 28 blood biomarkers revealed genes with heterogeneous effect sizes and directionality, highlighting the complex regulation of these traits [17].

• Rare Disease Diagnosis: The geneBurdenRD framework applied to 34,851 cases in the 100,000 Genomes Project identified 141 new disease-gene associations, five of which had independent supporting evidence, demonstrating the power of large-scale burden testing for rare Mendelian diseases [3].

This protocol provides a comprehensive framework for optimizing power in gene-based burden tests through rigorous variant quality control and strategic integration of functional annotations, enabling more powerful discovery of rare variant associations in complex and Mendelian diseases.

The identification of rare genetic variants associated with human diseases and traits represents a significant challenge in genetic association studies due to their low minor allele frequencies [18]. Single-variant tests often lack the statistical power to detect these associations, leading to the development of gene or region-based association tests that aggregate the effects of multiple rare variants [18]. As large-scale sequencing studies and biobanks continue to grow, meta-analysis—the combination of summary statistics across multiple cohorts—has emerged as a powerful strategy for enhancing the detection power of rare variant association tests [18] [41].

However, existing meta-analysis methods face substantial computational and statistical challenges. Methods such as MetaSTAAR can exhibit notably inflated type I error rates under imbalanced case-control ratios, a common scenario in biobank-based disease phenotype studies [18]. Additionally, the computational burden of constructing separate linkage disequilibrium (LD) matrices for each phenotype poses significant limitations for phenome-wide analyses [18]. To address these challenges, two novel methods have been developed: Meta-SAIGE and REMETA. These approaches offer improved type I error control and enhanced computational efficiency through innovative statistical and computational frameworks [18] [41].

This article explores the methodologies and applications of Meta-SAIGE and REMETA in the context of gene-based burden tests for rare variants, providing detailed protocols for implementation and comparing their performance against existing approaches.

Meta-SAIGE: Scalable and Accurate Rare Variant Meta-Analysis

Meta-SAIGE extends the SAIGE-GENE+ framework to meta-analysis, addressing key limitations of existing methods through several innovative approaches [18] [52]. The method employs two-level saddlepoint approximation (SPA) to control type I error inflation in the presence of case-control imbalance [18]. This includes SPA on score statistics of each cohort and a genotype-count-based SPA for combined score statistics from multiple cohorts [18]. Additionally, Meta-SAIGE allows the reuse of a single sparse LD matrix across all phenotypes, significantly reducing computational costs in phenome-wide analyses [18].

The Meta-SAIGE workflow consists of three main steps. First, per-variant level association summaries and a sparse LD matrix are prepared for each cohort using SAIGE, which derives score statistics (S) for both continuous and binary traits [18]. Second, score statistics from multiple cohorts are consolidated into a single statistic, with the covariance matrix of score statistics calculated in a sandwich form: Cov(S) = V^{1/2} Cor(G) V^{1/2}, where Cor(G) is the correlation matrix of genetic variants calculated from the sparse LD matrix Ω, and V is the diagonal matrix of the variance of S [18]. Finally, Meta-SAIGE performs rare variant tests—including Burden, SKAT, and SKAT-O set-based tests—utilizing various functional annotations and maximum minor allele frequency (MAF) cutoffs [18].

REMETA addresses the challenge of storing and sharing covariance matrices between variants for each study and trait of interest by introducing a reference LD matrix approach [41]. This method uses a single sparse covariance reference file per study that is rescaled for each phenotype using single-variant summary statistics, eliminating the need for separate LD matrices for each trait [41].

The REMETA workflow integrates seamlessly with the REGENIE software and involves three key steps [41]. First, reference LD matrices are constructed for each study, a process that needs to be performed only once per study [41]. Second, single-variant association testing is conducted using REGENIE on array genotypes for each study and trait, accounting for relatedness, population structure, and polygenicity [41]. Finally, gene-based meta-analysis is performed using REMETA, which applies a scaling and transformation to the reference LD files based on the trait being meta-analyzed, and computes burden tests, optimal sequence kernel association tests (SKATO), and aggregated Cauchy association tests (ACATV) [41].

Table 1: Key Features of Meta-SAIGE and REMETA

Feature	Meta-SAIGE	REMETA
Primary Innovation	Two-level saddlepoint approximation	Reference LD matrix rescaled per phenotype
Type I Error Control	Excellent for binary traits with case-control imbalance	Well-calibrated for binary traits with imbalance
Computational Efficiency	Reuses LD matrix across phenotypes	Single reference LD file per study
Software Integration	Extends SAIGE-GENE+	Integrates with REGENIE
Storage Requirements	O(MFK + MKP) for P phenotypes	Substantially reduced vs. trait-specific LD matrices
Key Applications	Low-prevalence binary traits, phenome-wide analyses	Large-scale exome sequencing studies, diverse cohorts

Experimental Validation and Performance

Type I Error Control

Simulation studies using UK Biobank (UKB) whole-exome sequencing (WES) data of 160,000 White British participants have demonstrated Meta-SAIGE's effectiveness in controlling type I error rates for binary traits with low prevalence [18]. When analyzing binary traits with a disease prevalence of 1% and a sample size ratio of 1:1:1 across three cohorts, methods without appropriate adjustment exhibited severe type I error inflation—approximately 100 times higher than the nominal level at α = 2.5 × 10^{-6} [18]. The application of SPA to score statistics of each cohort reduced but did not eliminate this inflation, while Meta-SAIGE's additional genotype-count-based SPA demonstrated well-controlled type I error rates across various scenarios [18].

REMETA has also shown accurate error control across a range of traits in the UK Biobank (n = 469,376), including body mass index (BMI), low-density lipoprotein (LDL), and various cancers with case-control ratios as extreme as 1:630 [41]. The method's approximation using reference LD matrices computed from all participants consistently produced accurate P values when applied to subsets of participants and different traits [41].

Statistical Power

Meta-SAIGE achieves statistical power comparable to joint analysis of individual-level data using SAIGE-GENE+ across various scenarios with different effect sizes of rare variants [18]. In contrast, the weighted Fisher's method, which aggregates SAIGE-GENE+ P values weighted by sample size, yields significantly lower power [18]. This demonstrates Meta-SAIGE's ability to enhance the detection of rare variant associations without requiring access to individual-level data from all cohorts.

In a practical application, Meta-SAIGE was used to conduct a meta-analysis of UK Biobank and All of Us WES data for 83 low-prevalence disease phenotypes, identifying 237 gene-trait associations at the exome-wide significance level [18] [52]. Notably, 80 of these associations were not significant in either dataset alone, underscoring the increased power afforded by this meta-analysis approach [18].

Computational Efficiency

Meta-SAIGE offers substantial computational advantages over existing methods. By reusing a single sparse LD matrix across all phenotypes, Meta-SAIGE significantly reduces storage requirements compared to approaches like MetaSTAAR that require separate LD matrices for each phenotype [18]. For meta-analyses of P different phenotypes across K cohorts with M variants, Meta-SAIGE requires O(MFK + MKP) storage, while MetaSTAAR requires O(MFKP + MKP) storage [18].

REMETA's approach of using precalculated reference LD matrices that can be shared across research groups further simplifies the computational landscape for large-scale meta-analyses [41]. The method employs a compact per-chromosome binary file format for efficiently storing and sharing LD matrices, which are indexed to allow fast access to the LD information of any gene [41].

Table 2: Performance Comparison of Meta-Analysis Methods

Performance Metric	Meta-SAIGE	REMETA	MetaSTAAR	Weighted Fisher's Method
Type I Error Control (Binary Traits)	Excellent	Well-calibrated	Inflated for imbalanced case-control ratios	Varies by implementation
Power	Comparable to individual-level analysis	Comparable to exact LD matrix	Comparable to individual-level analysis	Substantially lower
Computational Efficiency	High (reusable LD matrices)	High (reference LD matrices)	Lower (trait-specific LD matrices)	High (only P-values required)
Storage Requirements	Moderate	Reduced	High	Low
Ease of Implementation	Requires summary statistics and LD matrices	Integrates with REGENIE workflow	Requires summary statistics and trait-specific LD matrices	Simple (only P-values required)

Application Notes and Protocols

Protocol for Implementing Meta-SAIGE

Step 1: Preparation of Summary Statistics and LD Matrices

• For each cohort, use SAIGE to derive per-variant score statistics (S) for both continuous and binary traits [18].
• Compute the variance and association P values, applying SPA and efficient resampling for binary phenotypes based on minor allele count (MAC) [18].
• Generate a sparse LD matrix (Ω) representing the pairwise cross-product of dosages across genetic variants in the region to test [18].
• Store the summary statistics and LD matrices for each cohort in the required format for meta-analysis.

Step 2: Combination of Summary Statistics

• Combine score statistics from multiple cohorts into a single superset [18].
• For binary traits, recalculate the variance of each score statistic by inverting the P value obtained from SAIGE [18].
• Apply the genotype-count-based SPA to improve type I error control in the meta-analysis setting [18].
• Calculate the covariance matrix of score statistics using the formula: Cov(S) = V^{1/2} Cor(G) V^{1/2}, where Cor(G) is derived from the sparse LD matrix Ω, and V is the diagonal matrix of the variance of S [18].

Step 3: Gene-Based Association Testing

• Conduct Burden, SKAT, and SKAT-O set-based tests using the combined per-variant statistics and covariance matrix [18].
• Utilize various functional annotations and maximum MAF cutoffs to enhance detection power [18].
• Identify and collapse ultrarare variants (MAC < 10) to improve type I error control and power while reducing computational costs [18].
• Apply the Cauchy combination method to combine P values corresponding to different functional annotations and MAF cutoffs for each gene or region [18].

Protocol for Implementing REMETA

Step 1: LD Matrix Construction

• Construct reference LD matrices for each study using REMETA [41].
• This step needs to be performed only once for each study, regardless of the number of traits to be analyzed [41].
• Use the custom binary file format for efficient storage and sharing of LD matrices [41].
• For conditional gene-based tests, include both WES data and imputed variants across the genome in the LD matrix construction [41].

Step 2: Single-Variant Association Testing

• Run REGENIE Step 1 on array genotypes for each study and trait with appropriate covariates to account for relatedness, population structure, and polygenicity [41].
• Use the polygenic scores from REGENIE Step 1 as additional covariates in REGENIE Step 2 [41].
• Perform association testing of individual variants in the WES dataset for each phenotype using REGENIE Step 2 with the --htp flag to produce detailed summary statistic output required by REMETA [41].
• Include all polymorphic variants in this analysis without applying any filters on minor allele count to ensure comprehensive variant coverage for downstream gene-based tests [41].

Step 3: Gene-Based Meta-Analysis

• Use REMETA to perform gene-based meta-analysis with inputs including REGENIE summary statistic files for each trait and study, REMETA LD files for each study, gene set and variant annotation files, and an optional list of variants to condition on [41].
• REMETA will apply scaling and transformation to the reference LD files based on the trait being meta-analyzed [41].
• Compute burden tests, SKATO variance component tests, and ACATV using allele frequency bins specified by the user [41].
• Combine tests into an overall 'GENE P' P value for each gene to facilitate interpretation and prioritization [41].

Workflow for Meta-SAIGE and REMETA Implementation

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Research Reagents and Software Tools

Tool/Reagent	Function	Application Context
SAIGE Software	Per-variant association testing accounting for case-control imbalance and sample relatedness	Generation of input summary statistics for Meta-SAIGE [18]
REGENIE Software	Whole-genome regression for complex trait association analysis	Input generation for REMETA workflow [41]
Sparse LD Matrices	Storage of linkage disequilibrium information between genetic variants	Reference data for both Meta-SAIGE and REMETA [18] [41]
Variant Annotation Resources	Functional prediction of variant effects on protein structure	Grouping protein-damaging variants in gene-based tests [41]
Cauchy Combination Method	Combining P values from different functional annotations and MAF cutoffs	Meta-SAIGE gene-based test integration [18]
Binary File Format for LD Matrices	Efficient storage and sharing of LD information	REMETA implementation to reduce computational burden [41]

Meta-SAIGE and REMETA represent significant advances in rare variant meta-analysis, addressing critical limitations in type I error control and computational efficiency that have hampered previous approaches. Through innovative statistical methods such as two-level saddlepoint approximation and reference LD matrices, these methods enable robust, large-scale meta-analyses of rare variants across multiple cohorts and phenotypes.

The application of these methods to real-world datasets has demonstrated their practical utility, identifying novel gene-trait associations that were not detectable in individual cohorts. As biobanks and sequencing consortia continue to grow, Meta-SAIGE and REMETA will play increasingly important roles in unlocking the genetic architecture of complex diseases and traits, ultimately advancing drug target discovery and precision medicine initiatives.

Researchers implementing these methods should carefully consider their specific study designs, particularly for binary traits with extreme case-control imbalances, and follow the detailed protocols provided to ensure accurate and reproducible results. The continued development and refinement of meta-analysis methods will further enhance our ability to detect and interpret rare variant associations in diverse populations.

Gene-based burden tests are a fundamental tool in the analysis of exome-wide association studies (EWAS), where they function by aggregating multiple rare genetic variants within a gene into a single composite burden score. This approach increases statistical power for detecting associations between genes and complex traits or diseases, particularly when individual rare variants have insufficient effect sizes to be detected alone [53] [54]. A significant challenge in this paradigm is that researchers typically calculate and test a spectrum of burden scores across diverse annotation classes and frequency bins. This process generates multiple, highly correlated statistical tests, complicating both multiple testing correction and the biological interpretation of results [53] [54].

The Sparse Burden Association Test (SBAT) represents a methodological advancement designed to address this specific limitation. SBAT is a gene-based test that performs a joint analysis of the entire set of correlated burden scores simultaneously. The method operates under a foundational biological assumption: causal burden scores likely influence the trait in the same effect direction. By testing the joint set of scores, SBAT simultaneously assesses the overall significance of the model fit and identifies the subset of burden scores that most parsimoniously explain the observed association, thereby enhancing interpretability and power in specific scenarios [53] [54].

Table 1: Core Problem and SBAT Solution

Aspect	Standard Burden Test Challenges	SBAT Solution
Testing Framework	Multiple, separate tests for different burden scores	Joint testing of all burden scores concurrently
Correlation Handling	Correlated tests complicate multiple-testing correction	Built-in handling of correlation structure
Result Interpretation	Difficult to pinpoint the driving annotation	Selects the most explanatory set of burden scores
Genetic Assumption	Varies by test	Causal scores act in the same direction

SBAT Application Notes

Key Algorithmic Features and Workflow

The SBAT methodology incorporates several distinctive features. It employs non-negative least squares estimation, a computational approach that inherently respects the assumption of consistent effect direction for causal burden scores. This procedure yields a sparse solution, meaning it automatically selects a limited set of burden scores that contribute meaningfully to the model, reducing noise and focusing on the most biologically plausible signals. Furthermore, the SBAT framework integrates the processes of statistical testing and model selection, resolving the problem of having to run and correct for numerous separate hypothesis tests [53] [54].

The logical workflow of the SBAT method, from data preparation to result interpretation, is outlined in the diagram below.

Performance and Calibration

The SBAT method has undergone rigorous evaluation using simulated data, which confirmed that the test is well-calibrated. This means that under the null hypothesis (no true genetic association), the p-values it produces are uniformly distributed, protecting against false positives. More importantly, these simulations highlighted specific scenarios where SBAT demonstrates superior statistical power compared to existing gene-based tests, particularly when multiple correlated burden scores within a gene contain weak but consistent signals of association [53] [54].

The practical utility of SBAT has been demonstrated through large-scale application to 73 quantitative traits using exome sequencing data from the UK Biobank. This real-world analysis confirmed that SBAT serves as a valuable addition to the genomicist's toolkit. It is most effectively used not as a standalone test, but as a complementary method alongside other established gene-based tests, as it can detect associations that might be missed by other approaches [53] [54].

Table 2: SBAT Performance and Application Context

Evaluation Aspect	Key Finding	Implication for Research
Type I Error Control	Test is well-calibrated in simulations [53]	Robust against false positives
Statistical Power	Outperforms existing tests in specific scenarios [53]	Increased detection rate for certain genetic architectures
Real-World Application	Successfully applied to 73 UK Biobank traits [53] [54]	Validated for use with large-scale biobank data
Recommended Use	A valuable additional test [53] [54]	Use complementarily with other gene-based tests

Experimental Protocols

Software Implementation and Data Preparation

The SBAT method is implemented within the widely used REGENIE software, making it accessible to the broader research community without the need for custom programming [53] [54]. Researchers must first ensure they have installed a recent version of REGENIE that includes the SBAT functionality.

Input data preparation requires two primary components. First, genotype data must be processed, typically from exome or whole-genome sequencing, and formatted into a .vcf or .bgen file. This data should be restricted to rare variants (e.g., minor allele frequency < 0.01) within the gene(s) of interest. Second, a phenotype file is needed, containing the trait of interest for all samples, along with any relevant covariates (e.g., age, sex, principal components for ancestry). The following code block illustrates a typical REGENIE command for running the SBAT analysis:

Step-by-Step Analysis Procedure

The analytical procedure for SBAT can be broken down into distinct stages, as visualized in the following workflow.

• Data Preprocessing: Begin by generating a set of burden scores for each gene. This involves grouping rare variants (e.g., MAF < 1%) by gene and creating multiple scores based on different functional annotations (e.g., loss-of-function, missense) and frequency bins [53] [54].
• Command Execution: Execute the REGENIE command with the --sbat flag, as shown in the code block above. The software will automatically handle the joint testing of all specified burden scores for each gene.
• Output Examination: The primary output includes a list of genes with a single, overall p-value from the joint test. Review the results to identify genes that surpass a multiple-testing-corrected significance threshold (e.g., exome-wide significance).
• Post-Hoc Interpretation: For significant genes, inspect the specific burden scores that were selected by the non-negative least squares algorithm. This sparse set represents the variant classes most likely driving the association.
• Downstream Validation: Significant findings from SBAT should be followed by replication in an independent cohort and investigated through functional genomic experiments to elucidate the biological mechanism.

The Scientist's Toolkit: Research Reagent Solutions

The following table details key resources and information required to implement an SBAT analysis effectively.

Table 3: Essential Materials and Resources for SBAT Analysis

Item / Resource	Function / Description	Implementation Note
REGENIE Software	The software platform implementing the SBAT test.	Ensure version compatibility and include the --sbat flag in the command.
Exome/Genome Sequencing Data	Source of rare genetic variants for burden score calculation.	Data should be in .vcf or .bgen format; quality control (QC) is critical.
Phenotype & Covariate Files	Contains trait measurements and adjustment variables.	Properly format according to REGENIE requirements; handle missing data.
Variant Annotation Tool (e.g., VEP)	Classifies variants by function (LOF, missense, etc.).	Used in pre-processing to define the annotation classes for burden scores.
High-Performance Computing (HPC) Cluster	Computational environment for running analyses.	REGENIE-SBAT on genome-wide scale requires substantial memory and compute time.
Burden Score Definitions	The specific rules for aggregating variants into a single score.	Pre-define multiple scores (e.g., based on frequency and predicted consequence).

Benchmarking Performance and Translating Insights into Clinical Value

Gene-based burden tests have become a cornerstone in rare variant association analysis, offering a powerful strategy to uncover the genetic underpinnings of complex diseases by aggregating signals from multiple rare variants within a gene or genomic region [41] [21]. The fundamental premise of these tests is to compensate for the low statistical power of single-variant tests when analyzing rare genetic variants with minor allele frequencies (MAFs) below 1% [21]. By combining evidence across multiple variants, burden tests can significantly enhance the discovery of novel gene-trait associations, particularly in the context of large-scale exome-wide association studies (ExWAS) and drug target discovery [41].

As the field of genetic association studies evolves with ever-increasing sample sizes from biobanks and consortia, the comparative performance of different gene-based tests—specifically their statistical power to detect true associations and their ability to control type I error rates (false positives)—has emerged as a critical consideration for researchers [55] [56]. The performance characteristics of these methods are not uniform; they vary substantially depending on the underlying genetic architecture, sample size, and specific methodological approaches [21]. This protocol provides a systematic framework for evaluating leading gene-based burden tests, with particular emphasis on their operating characteristics in diverse research scenarios relevant to scientists and drug development professionals.

Comparative Power Analysis of Gene-Based Burden Tests

The statistical power of gene-based burden tests is influenced by multiple factors, including the proportion of causal variants within a gene, their effect sizes and directions, the linkage disequilibrium (LD) pattern between variants, and the total sample size [55] [21]. Understanding how these factors interact with different methodological approaches is essential for selecting the optimal test in a given research context.

Table 1: Comparative Power of Gene-Based Tests Under Different Genetic Architectures

Test Method	High Proportion of Causal Variants	Mixed Effect Directions	Few Causal Variants	Key Assumptions
Burden Tests	High power [21]	Reduced power [26]	Moderate power [26]	All causal variants have same effect direction [26]
Variance Component Tests (SKAT)	Moderate power [26]	High power [26]	Good power [57]	Variants can have different effect directions [26]
Combined Tests (SKAT-O)	High power [18]	High power [18]	Good power [18]	Adaptive combination of burden and variance components [18]
Omnibus Tests (ACAT-V)	Good power [26]	Good power [26]	High power [26]	Combines multiple test statistics; robust to different architectures [26]

Recent methodological advances have introduced sophisticated approaches that jointly model multiple burden scores. The Sparse Burden Association Test (SBAT) simultaneously assesses multiple burden scores across different frequency thresholds and annotation classes under the assumption that causal burden scores act in the same effect direction [26]. This method provides both model inference with controlled type I error and model selection, identifying which frequency bins and annotation classes harbor causal variants [26]. Similarly, the Overall method aggregates information from burden tests, SKAT, and SKAT-O with multiple eQTL-derived weights using an extended Simes procedure, demonstrating enhanced power across diverse genetic architectures [57].

For binary traits with case-control imbalance, methods like Meta-SAIGE employ saddlepoint approximation (SPA) to maintain power while controlling type I error [18]. This is particularly relevant for biobank-based studies of low-prevalence diseases, where standard methods may suffer from substantial power loss [18].

Table 2: Power Determinants in Rare Variant Association Studies

Factor	Impact on Power	Optimal Method
Proportion of Causal Variants	Aggregation tests more powerful only when substantial proportion are causal (>55% in some scenarios) [21]	Burden tests when high proportion; single-variant or ACAT-V when low proportion [21]
Effect Direction	Burden tests powerful when all effects in same direction; variance component tests when mixed directions [26]	SKAT or SKAT-O for mixed directions; burden tests for uniform directions [26]
Sample Size	Aggregation tests require large samples (n=100,000+ for 0.1% heritability) [21]	Meta-analysis methods (e.g., REMETA, Meta-SAIGE) for combining cohorts [41] [18]
Case-Control Imbalance	Severe imbalance reduces power for binary traits [18]	Saddlepoint approximation methods (e.g., SAIGE, Meta-SAIGE) [18]
Variant Annotation	Functional annotation improves power by focusing on deleterious variants [21]	Methods incorporating functional weights (e.g., SBAT, Overall) [26] [57]

Type I Error Control in Gene-Based Tests

Proper control of type I error rates (false positives) is equally important as maintaining power when comparing gene-based tests. Different methods exhibit varying performance in preserving nominal error rates under challenging conditions such as high correlation between variants (linkage disequilibrium), case-control imbalance, and inclusion of rare variants.

Brown's method has been shown to be among the most robust techniques for controlling type I error rates when combining p-values from correlated variants within a gene [55]. Simulation studies based on HapMap project data demonstrated that this method effectively maintains error rates at the desired nominal level across genes with different characteristics and LD patterns [55]. In contrast, Fisher's method, which assumes independence between variants, typically leads to inflated type I error rates when applied to genetic variants in LD [55].

For rare variant analysis in binary traits with substantial case-control imbalance, methods like Meta-SAIGE that employ two-level saddlepoint approximation (SPA) have demonstrated superior type I error control compared to alternatives [18]. In simulations with disease prevalence of 1% and sample size ratio of 1:1:1 across three cohorts, methods without appropriate adjustment showed type I error rates nearly 100 times higher than the nominal level at α = 2.5×10⁻⁶, whereas Meta-SAIGE with GC-based SPA maintained well-controlled error rates [18].

The GATES method provides another approach for maintaining proper error control through an extended Simes procedure that accounts for the effective number of independent SNPs in a gene, estimated from the eigenvalues of the LD matrix [58] [59]. This method does not require permutation or simulation to evaluate empirical significance, making it computationally efficient while controlling type I error rates regardless of gene size and LD pattern [59].

Multivariate gene-based tests that analyze multiple correlated phenotypes simultaneously present additional challenges for error control. Studies have shown that combinations of multivariate association methods with gene-based tests can lead to inflated type I errors if not properly calibrated. For instance, TATES and MultiPhen paired with VEGAS demonstrate inflated type I error across various scenarios, while these same multivariate methods paired with GATES maintain correct error rates [58].

Experimental Protocols

Protocol 1: Gene-Based Association Analysis Using REGENIE/REMETA Workflow

The REGENIE/REMETA pipeline provides a computationally efficient approach for gene-based meta-analysis of rare variants using summary statistics [41]. This protocol details the three-step workflow for analyzing multiple traits across large-scale genetic studies.

Step 1: LD Matrix Construction

• Construct reference LD matrices for each study using the REMETA software [41]
• Generate sparse covariance matrices from the whole-exome sequencing (WES) dataset
• Apply the QR decomposition to retain only linearly independent columns for numerical stability [26]
• Store matrices in compact per-chromosome binary file format for efficient storage and retrieval [41]

Step 2: Single-Variant Association Testing

• Run REGENIE Step 1 on array genotypes for each study and trait with appropriate covariates to account for relatedness, population structure, and polygenicity [41]
• Use the resulting polygenic scores as additional covariates in REGENIE Step 2
• Perform association testing of individual variants in the WES dataset for each phenotype
• Apply the --htp flag in REGENIE Step 2 to produce detailed summary statistic output required by REMETA [41]
• Critical: Include all polymorphic variants without minor allele count filters to enable comprehensive downstream gene-based tests [41]

Step 3: Meta-Analysis

• Input REGENIE summary statistics, REMETA LD files, gene sets, and variant annotation files into REMETA [41]
• Specify an optional list of variants to condition analysis upon
• Apply scaling and transformation to reference LD files based on the specific trait being analyzed [41]
• Compute burden tests, SKATO variance component tests, and ACATV using user-specified allele frequency bins [41]
• Generate overall 'GENE P' values for each gene by combining evidence across different tests [41]

Workflow for REGENIE/REMETA Gene-Based Analysis

Protocol 2: Type I Error Assessment for Binary Traits with Case-Control Imbalance

This protocol outlines the procedure for evaluating type I error control in gene-based tests, particularly for challenging scenarios with highly imbalanced case-control ratios, as commonly encountered in biobank studies of rare diseases.

Simulation Design:

• Utilize real whole-exome sequencing data from large biobanks (e.g., UK Biobank) as genotype foundation [18]
• Generate null phenotypes with varying prevalences (e.g., 5%, 1%) under the null hypothesis of no association [18]
• Divide samples into multiple non-overlapping cohorts with different size ratios (e.g., 1:1:1, 4:3:2) to mimic multi-cohort meta-analysis [18]
• Conduct a minimum of 60 simulation replicates with approximately 1 million tests per replicate to achieve sufficient precision for evaluating exome-wide significance thresholds [18]

Method Application:

• Apply target gene-based tests (e.g., Meta-SAIGE, MetaSTAAR) to simulated data [18]
• For methods requiring type I error correction (e.g., Meta-SAIGE), implement saddlepoint approximation (SPA) on score statistics of each cohort [18]
• Apply genotype-count-based SPA for combined score statistics from multiple cohorts [18]
• Compare unadjusted approaches to those with SPA and GC-based SPA adjustments [18]

Error Rate Calculation:

• Calculate empirical type I error rates as the proportion of tests rejecting the null hypothesis at specified significance thresholds (e.g., α = 2.5×10⁻⁶) [18]
• Compare empirical rates to nominal levels to assess potential inflation or deflation
• Evaluate consistency of error control across different prevalence scenarios and cohort configurations [18]

Protocol 3: Power Comparison Between Single-Variant and Aggregation Tests

This protocol describes a framework for determining when aggregation tests provide power advantages over single-variant tests for rare variant association analysis.

Analytical Power Calculations:

• Use available online tools (e.g., R Shiny app: https://debrajbose.shinyapps.io/analytic_calculations/) for initial power estimations [21]
• Assume an additive genetic model with c causal out of v total rare variants in a gene [21]
• Specify region heritability (h²), sample size (n), and variant characteristics [21]
• Calculate non-centrality parameters for single-variant, burden test, and SKAT statistics under the assumption of independent variants [21]

Simulation-Based Power Assessment:

• Simulate traits based on real genotype data from large-scale sequencing studies (e.g., UK Biobank with 378,215 unrelated participants) [21]
• Vary key parameters: proportion of causal variants (5%, 15%, 55%), effect sizes, and effect directions [21]
• Consider different variant masks focusing on likely high-impact variants (e.g., protein-truncating variants, deleterious missense variants) [21]
• For each scenario, apply both single-variant and aggregation tests to the same simulated data [21]

Power Comparison:

• Compute power as the proportion of simulations where each test rejects the null hypothesis at the specified significance level
• Identify scenarios where aggregation tests outperform single-variant tests and vice versa
• Determine thresholds for the proportion of causal variants where aggregation tests become advantageous [21]

The Scientist's Toolkit

Table 3: Essential Research Reagents and Computational Tools for Gene-Based Burden Tests

Tool/Resource	Function	Application Context
REGENIE/REMETA [41] [26]	Efficient gene-based testing and meta-analysis	Large-scale exome-wide association studies; meta-analysis across diverse cohorts
Meta-SAIGE [18]	Rare variant meta-analysis with accurate type I error control	Binary traits with case-control imbalance; phenome-wide analyses
SBAT [26]	Joint testing of multiple burden scores with sparsity	Quantitative traits; identifying relevant frequency bins and annotation classes
GATES [58] [59]	Rapid gene-based testing using extended Simes procedure	Genome-wide association studies with common variants; efficient p-value combination
Overall Method [57]	Omnibus test combining multiple gene-based tests with eQTL weights	Integrating functional genomic data; enhancing power through multi-test aggregation
1000 Genomes Project [58]	Reference panel for LD estimation	LD estimation when individual-level data unavailable
GTEx Portal [57]	Source of eQTL-derived weights	Incorporating functional information into gene-based tests
UK Biobank WES Data [21] [18]	Large-scale reference dataset	Method validation; simulation studies; power calculations

The comparative analysis of power and type I error rates presented in this protocol reveals that method selection for gene-based burden tests should be guided by the specific research context, particularly the underlying genetic architecture and study design. Burden tests demonstrate optimal performance when a substantial proportion of aggregated variants are causal and act in the same effect direction, while variance component tests like SKAT maintain advantage when causal variants have mixed effect directions [21] [26]. Omnibus tests such as ACAT-V and Overall provide robust performance across diverse genetic architectures by combining evidence from multiple testing approaches [26] [57].

For binary traits with case-control imbalance, methods incorporating saddlepoint approximation (e.g., Meta-SAIGE) are essential for maintaining proper type I error control without sacrificing power [18]. In large-scale meta-analyses across multiple cohorts, computationally efficient approaches like REMETA that reuse reference LD matrices across phenotypes offer substantial practical advantages while maintaining statistical accuracy [41].

As rare variant association studies continue to scale with expanding biobanks and consortia, the thoughtful selection of gene-based burden tests informed by their power and error characteristics will remain crucial for advancing our understanding of complex disease genetics and identifying novel therapeutic targets.

The field of human genetics has long relied on two distinct approaches for uncovering the genetic basis of traits and diseases: genome-wide association studies (GWAS) that identify common variants with small effects, and rare-variant burden tests that detect aggregated signals from rare variants with larger effects. While these methods are conceptually similar, they systematically prioritize different genes, raising fundamental questions about optimal gene prioritization strategies [60]. Historically, common and rare variant analyses were conducted independently, but emerging evidence demonstrates that integrating these approaches provides a more comprehensive understanding of disease architecture and enables superior risk prediction. This protocol examines the theoretical foundation, empirical validation, and practical implementation of integrated variant models for complex disease prediction.

Theoretical Foundation: Why Integration Works

Complementary Biological Insights from Different Variant Classes

Common and rare variants provide distinct yet complementary biological insights. GWAS prioritize genes near trait-specific variants, whereas burden tests prioritize trait-specific genes [60]. This difference arises because non-coding variants identified by GWAS can be context-specific, allowing them to identify highly pleiotropic genes, while burden tests generally cannot. Importantly, both study designs are affected by distinct trait-irrelevant factors, including gene length and random genetic drift, further complicating their interpretation [60].

The conceptual framework for gene prioritization involves two key criteria:

• Trait importance: How much a gene quantitatively affects a trait
• Trait specificity: The importance of a gene for the trait under study relative to its importance across all traits [60]

Burden tests naturally favor trait-specific genes, while GWAS captures both specific and pleiotropic genes, explaining their differential performance in gene discovery.

Power Considerations for Variant Detection

The relative power of aggregation tests versus single-variant tests depends critically on the underlying genetic architecture. Analytical calculations and simulations reveal that aggregation tests are more powerful than single-variant tests only when a substantial proportion of variants are causal, with power strongly dependent on the specific genetic model and set of rare variants aggregated [21]. For example, when aggregating protein-truncating variants (PTVs) and deleterious missense variants, aggregation tests become more powerful than single-variant tests for >55% of genes when PTVs, deleterious missense variants, and other missense variants have 80%, 50%, and 1% probabilities of being causal, respectively (with n=100,000 and h²=0.1%) [21].

Table 1: Conditions Favoring Aggregation Tests Over Single-Variant Tests

Factor	Condition Favoring Aggregation Tests	Typical Threshold
Proportion of causal variants	High proportion	>25-50% depending on effect sizes
Sample size	Large cohorts	n > 50,000
Region heritability	Moderate to high	h² > 0.1%
Variant impact	High-impact variants	PTVs, deleterious missense
Effect direction	Consistent direction	All risk-increasing or all protective

Empirical Evidence from Disease Applications

Breast Cancer Risk Stratification

A landmark study of 25,591 women in the Healthy Nevada Project demonstrated the superior performance of integrated risk assessment. The combined approach identified 11.4% of women as high-risk based on possessing either a predicted loss-of-function (pLOF) variant in BRCA1, BRCA2, or PALB2 (hazard ratio = 10.4); a pLOF variant in ATM or CHEK2 (hazard ratio = 3.4); or being in the top 10% of the polygenic risk score (PRS) distribution (hazard ratio = 2.4) [61].

Critically, the study revealed important interaction effects between rare variants and polygenic risk. Women with a pLOF in ATM or CHEK2 who also ranked in the top 50% of the PRS displayed dramatically elevated risk (39.2% probability of breast cancer at age 70), whereas their counterparts in the bottom 50% of the PRS were not at high risk (14.4% probability at age 70) [61]. This gene-environment interaction (GxE) equivalent at the genetic level highlights the necessity of integrated models.

Table 2: Performance of Genetic Screening Approaches for Breast Cancer

Screening Approach	High-Risk Group Size	Hazard Ratio (95% CI)	Key Findings
Family history only	Not reliably ascertained	N/A	78% of family history records occurred after diagnosis
Rare variants (BRCA1/2, PALB2)	~1-2%	10.4 (8.1-13.5)	High penetrance but limited population impact
Rare variants (ATM, CHEK2)	~2-3%	3.4 (2.4-4.8)	Moderate penetrance, highly dependent on PRS
Polygenic risk score (top 10%)	10%	2.4 (2.0-2.8)	Moderate population-wide impact
Combined model	11.4%	Varied by profile	Optimized sensitivity and specificity

Atrial Fibrillation Genomic Risk Prediction

A comprehensive cohort study using whole-genome sequence data from 416,085 UK Biobank participants developed an integrated genomic model for atrial fibrillation (IGM-AF) incorporating polygenic risk score (PRS), a composite rare variant gene set (AFgeneset), and somatic variants associated with clonal hematopoiesis of indeterminate potential (CHIP) [62].

Each component independently predicted incident AF:

• PRS: HR per 1 SD = 1.65 (95% CI: 1.63-1.67)
• AFgeneset: HR = 1.63 (95% CI: 1.52-1.75)
• CHIP: HR = 1.26 (95% CI: 1.15-1.38)

The 5-year cumulative incidence of AF was at least 2-fold higher among individuals with all three genetic drivers (common, rare, and somatic) compared to those with only one driver [62]. Integration of IGM-AF with a clinical risk model (CHARGE-AF) showed higher predictive performance (C statistic = 0.80) compared to IGM-AF and CHARGE-AF alone, with improved reclassification (net reclassification index = 0.08) [62].

Enhanced Rare Disease Diagnosis with AI

The popEVE AI model represents a breakthrough in rare variant interpretation, combining deep evolutionary information with human population data to predict variant pathogenicity [63]. Unlike previous models that could not easily compare variants across different genes, popEVE produces a calibrated score for each variant that can be compared across genes, enabling prioritization of clinically relevant findings.

In validation studies involving approximately 30,000 patients with severe developmental disorders who had not received diagnoses, popEVE led to diagnoses in about one-third of cases and identified variants in 123 genes not previously linked to developmental disorders [63]. The model demonstrated 15-fold enrichment for highly deleterious variants in affected individuals compared with controls, with 98% accuracy in ranking likely causal variants above all inherited variants in the same patient [64].

Integrated Analysis Protocol

Meta-Analysis of Rare Variants Across Cohorts

For large-scale rare variant association testing, Meta-SAIGE provides a scalable method for rare variant meta-analysis that accurately controls type I error rates, particularly for low-prevalence binary traits [18]. The protocol involves three principal steps:

Step 1: Cohort-Level Preparation

• Use SAIGE to derive per-variant score statistics (S) for both continuous and binary traits
• Calculate variance and association P values, adjusting for case-control imbalance and sample relatedness
• Generate a sparse linkage disequilibrium (LD) matrix (Ω) representing pairwise cross-products of dosages across genetic variants in the region of interest
• For binary phenotypes, compute P values using saddlepoint approximation (SPA) for variants with minor allele count (MAC) ≥ 20 and efficient resampling for MAC < 20 [18]

Step 2: Summary Statistics Combination

• Combine score statistics from multiple cohorts into a single superset
• For binary traits, recalculate the variance of each score statistic by inverting the P value from SAIGE
• Apply genotype-count-based SPA to improve type I error control
• Calculate the covariance matrix of score statistics using the formula: Cov(S) = V¹ᐧ²Cor(G)V¹ᐧ², where Cor(G) is the correlation matrix from the sparse LD matrix Ω, and V is the diagonal variance matrix [18]

Step 3: Gene-Based Association Testing

• Conduct Burden, SKAT, and SKAT-O set-based tests using various functional annotations and MAF cutoffs
• Identify ultrarare variants (MAC < 10) and collapse them to enhance power while controlling type I error
• Use the Cauchy combination method to combine P values corresponding to different functional annotations and MAF cutoffs for each gene or region [18]

Integrated Risk Model Development

Stage 1: Component Generation

• Polygenic Risk Score: Develop PRS using pruning and thresholding, LDpred, or Bayesian polygenic modeling based on large-scale GWAS summary statistics
• Rate Variant Burden: Aggregate qualifying rare variants (typically PTVs and deleterious missense with MAF < 0.1%) within established disease-associated genes using validated pathogenicity prediction tools
• Somatic Variant Assessment: Identify clonal hematopoiesis (CHIP) variants from blood-derived DNA sequencing data using variant allele fraction thresholds (typically >2%)

Stage 2: Model Integration

• Apply multivariable Cox proportional hazards or logistic regression models with clinical covariates (age, sex, principal components)
• Include PRS (as continuous Z-score), rare variant carrier status (binary), and somatic variant status (binary) as predictors
• Test for interactions between PRS and rare variant carrier status using multiplicative interaction terms

Stage 3: Validation and Calibration

• Perform internal validation using bootstrapping or cross-validation
• Conduct external validation in independent cohorts when available
• Assess calibration using observed vs. expected event rates across risk deciles
• Calculate discrimination statistics (C-index, AUC) and reclassification metrics (NRI, IDI)

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Resources for Integrated Variant Studies

Resource Category	Specific Tools/Methods	Application Note
Variant Calling	GATK, DeepVariant, GLnexus	Joint calling improves rare variant detection; consider cohort-specific batch effects
Pathogenicity Prediction	popEVE, EVE, REVEL, AlphaMissense	popEVE combines evolutionary and population constraint for cross-gene comparison [63]
Burden Test Implementation	SAIGE-GENE+, Meta-SAIGE, STAAR	Meta-SAIGE controls type I error in meta-analysis of binary traits [18]
Rare Variant Association Tests	Burden, SKAT, SKAT-O, ACAT	Choose based on genetic architecture; SKAT-O optimal for mixed effect directions [65] [66]
Polygenic Risk Scores	PRSice, LDpred, lassosum, CT-SLEB	Ensure ancestry-matched effect size estimates; consider portability across populations
Meta-Analysis Tools	Meta-SAIGE, RAREMETAL, METAL	Meta-SAIGE reuses LD matrices across phenotypes for computational efficiency [18]
Functional Annotation	ANNOVAR, VEP, FunSeq2	Integrate regulatory annotations for non-coding variants in GWAS loci
Population Frequency Databases	gnomAD, UK Biobank, All of Us	Use ancestry-specific frequencies for rare variant filtering

Discussion and Future Directions

The integration of common and rare variants represents a paradigm shift in complex disease genetics, moving beyond the traditional dichotomy that treated these variant classes separately. Evidence across multiple diseases consistently demonstrates that integrated models provide superior risk prediction compared to approaches using either variant class alone [61] [62]. This improved performance stems from the biological complementarity of common and rare variants, with common variants capturing distributed polygenic risk and rare variants identifying high-effect disruptions in biologically central pathways.

Methodological advances continue to address key challenges in variant integration. For rare variant analysis, methods like Meta-SAIGE overcome the historical limitations of type I error inflation in meta-analyses of binary traits [18]. Artificial intelligence approaches like popEVE enhance rare variant interpretation by combining evolutionary and population genetic data to enable cross-gene comparison of variant deleteriousness [63]. For common variants, increasingly sophisticated PRS methods improve cross-ancestry portability and predictive accuracy.

Future directions in the field include:

• Development of ancestry-aware integrated models to ensure equitable translation across diverse populations
• Incorporation of non-coding regulatory variants through transcriptome-based methods
• Integration of somatic variation in non-blood tissues for broader disease applications
• Longitudinal modeling of genetic risk across the lifespan
• Clinical trial frameworks to evaluate interventions based on integrated genetic risk

The protocol outlined herein provides a framework for researchers to implement integrated variant analyses, with applications spanning disease gene discovery, risk prediction, and therapeutic target identification. As genetic datasets continue to grow in size and diversity, the integration of common, rare, and somatic variants will increasingly become the standard approach for unraveling the genetic architecture of complex diseases and translating these insights into clinical practice.

Gene-based burden tests represent a cornerstone of modern rare variant analysis, providing a powerful framework for detecting associations between aggregated rare genetic variants and complex phenotypes. The advent of large-scale biobanks, notably the UK Biobank (UKBB) and the All of Us (AoU) Research Program, has revolutionized the scale and scope of these analyses by providing unprecedented sample sizes and diverse genetic data. These resources enable researchers to move beyond traditional common variant analyses to explore the substantial contribution of rare coding and noncoding variants to human disease. The integration of whole-genome sequencing (WGS) data with deep phenotypic information allows for systematic investigation of rare variants across the allele frequency spectrum, revealing new biological mechanisms and potential therapeutic targets [67] [36]. This application note details standardized protocols for conducting gene-based burden analyses within these biobanks, highlighting key applications, methodological considerations, and validation frameworks essential for robust genetic discovery.

Comparative Biobank Characteristics

Table 1: Design and Genomic Characteristics of Major Biobanks

Biobank Feature	UK Biobank (UKBB)	All of Us (AoU)
Total Participants	~500,000	245,388 (WGS released, goal: >1 million)
WGS Data Available	490,640 participants	245,388 participants (as of Feb 2024)
Ancestry Diversity	93.5% European, 6.5% other (African, South Asian, East Asian)	77% from groups underrepresented in biomedical research
Key Phenotypic Data	EHRs, surveys, physical measures, imaging	EHRs, surveys, physical measures, wearable data
Exome Data	469,818 participants	Integrated with WGS data
Notable Features	Extensive phenotyping, longitudinal follow-up	Focus on health equity, diverse recruitment

Quantitative Insights from Burden Analyses

Table 2: Representative Findings from Gene-Based Burden Analyses

Phenotype	Key Genes Identified	Effect Size/Direction	Biobank	Sample Size
Hyperlipidaemia	LDLR, PCSK9, APOC3, ANGPTL3, ABCG5, NPC1L1	Increased risk (LDLR, ABCG5); Protective (others)	UKBB	470,000 exomes
Body Mass Index	RIF1, UBR3, TNK2	PTVs associated with higher BMI (e.g., +2.66 kg/m² for RIF1)	UKBB + AoU	489,941 UKBB + ~220,000 AoU
Type 2 Diabetes	IRS2, UBR3, NAA15, RMC1	PTVs increase risk (e.g., OR=6.4 for IRS2)	UKBB + AoU	489,941 UKBB + ~220,000 AoU
Aortic Stenosis	LDLR, APOB, PCSK9	Rare variants associated with increased risk	AoU	CDRv7 dataset
Chronic Kidney Disease	APOL3 (in epistasis with APOL1)	PTV increases risk	AoU	CDRv7 dataset

Experimental Protocols for Gene-Based Burden Analysis

Variant Annotation and Quality Control

Purpose: To classify rare genetic variants based on functional impact and ensure data quality prior to association testing.

Procedure:

• Variant Annotation: Annotate all variants using Ensembl's Variant Effect Predictor (VEP) with additional functional prediction scores (REVEL, PolyPhen, SIFT) [68] [36].
• Variant Categorization: Classify variants into:
  • Protein-truncating variants (PTVs): Stop-gain, frameshift, essential splice-site
  • Damaging missense: REVEL score >0.5 or >0.7
  • Other missense: REVEL score ≤0.5
  • Synonymous: Used as negative controls [36]
• Frequency Filtering: Retain variants with minor allele frequency (MAF) <0.1% (for very rare) or <1% (for rare) in population controls [36].
• Quality Filters: Exclude variants with:
  • >10% missing genotypes in cases or controls
  • Deviation from Hardy-Weinberg equilibrium (p<1×10⁻⁶)
  • Poor sequencing quality metrics [68]

Burden Score Construction

Purpose: To aggregate rare variants within genes to increase statistical power for association testing.

Procedure:

• Variant Weighting Scheme: Assign weights based on:
  • Frequency: Parabolic weighting where MAF=0.01 gets weight=1 and MAF≈0 gets weight=10 [68]
  • Functional Impact:
    • PTVs: weight=100
    • Damaging missense: weight=25-30 (5 base + 20 for SIFT deleterious + 5 for PolyPhen damaging)
    • Other missense: weight=5 [68]
• Gene-Based Aggregation: For each sample, calculate burden score as:
  • ( \text{Burden Score}g = \sum{i=1}^{n} w{f,i} \times w{a,i} \times Ci )
  • Where ( w{f,i} ) = frequency weight, ( w{a,i} ) = functional annotation weight, and ( Ci ) = allele count for variant i [68]
• Mask Definition: Create variant masks for testing:
  • PTV-only mask
  • PTV + damaging missense (REVEL>0.7)
  • PTV + damaging missense (REVEL>0.5) [36]

Statistical Association Testing

Purpose: To test for significant associations between gene burden scores and phenotypes while controlling for confounding variables.

Procedure:

• Regression Framework: Implement mixed-effects or generalized linear models:
  • For continuous traits: Linear regression
  • For binary traits: Logistic regression [69] [68]
• Covariate Adjustment: Include as covariates:
  • Age, sex, age²
  • Genotyping platform/center
  • First 20 genetic principal components (for ancestry)
  • Relatedness (using genetic relationship matrix) [36] [70]
• Significance Thresholding: Apply exome-wide significance threshold:
  • ( P < 6.15 \times 10^{-7} ) for ~81,350 tests (19,457 genes × 4 masks) [36]
  • Bonferroni correction for hypothesis-driven analyses of specific gene sets
• Replication Framework: Test significant associations in independent datasets:
  • Cross-biobank validation (UKBB → AoU)
  • Ancestry-specific replication [69] [36]

Signaling Pathways and Analytical Workflows

Diagram 1: Gene-based burden analysis workflow showing key stages from data preparation through biological interpretation.

Diagram 2: IRS2 insulin signaling pathway showing mechanistic links to type 2 diabetes and chronic kidney disease risk.

The Scientist's Toolkit: Research Reagent Solutions

Table 3: Essential Analytical Tools and Resources for Burden Testing

Tool/Resource	Function	Application Example	Source/Reference
REGENIE with SBAT	Joint testing of multiple burden scores	Sparse Burden Association Test for correlated burden scores	[71]
SCOREASSOC	Weighted burden analysis	Implementing frequency and function-based variant weights	[68]
GENEVARASSOC	Variant annotation for burden testing	Generating input files for weighted burden analysis	[68]
Genebass Browser	Rapid exploration of rare-variant associations	Systematic scanning of 4,529 phenotypes in UKBB	[72]
Variant Effect Predictor (VEP)	Functional variant annotation	Categorizing PTVs, missense, and regulatory variants	[68] [70]
REVEL Scores	Missense variant deleteriousness	Defining damaging missense variants (threshold >0.5/0.7)	[36]
SF-Relate	Privacy-preserving genetic relatedness	Identifying genetic relatives across distributed datasets	[73]

Applications and Case Studies

Cross-Biobank Validation for Cardiometabolic Traits

The integration of UK Biobank and All of Us data has proven particularly powerful for validating novel genetic associations with cardiometabolic traits. Recent WGS analysis of 489,941 UKBB participants identified protein-truncating variants (PTVs) in IRS2 with a substantial effect on type 2 diabetes risk (odds ratio 6.4, 34% case prevalence among carriers). These findings were replicated in approximately 220,000 individuals from the All of Us study, confirming IRS2's role in diabetes pathogenesis [36]. Similarly, PTVs in UBR3 demonstrated independent effects on both body mass index and type 2 diabetes risk across both biobanks, highlighting pleiotropic effects that may inform therapeutic development.

Ancestry-Inclusive Discovery Efforts

The emphasis on diverse recruitment in AoU (77% from populations underrepresented in biomedical research) has enabled more inclusive genetic discovery. A pan-ancestry analysis of sequencing data from 748,879 individuals across multiple biobanks, including All of Us, identified 363 significant gene-disease associations, including novel associations for UBR3 in cardiometabolic disease and YLPM1 in psychiatric disease [69]. This study demonstrated that effect sizes for rare protein-disrupting variants were largely concordant across ancestries (βDeming = 0.7-1.0), supporting the transferability of genetic findings across populations when adequate representation is achieved.

Noncoding Variant Discovery Through WGS

The comprehensive nature of whole-genome sequencing in both UKBB and AoU enables investigation beyond coding regions. A recent analysis of WGS data from ~50,000 UKBB participants identified 604 independent rare noncoding single-variant associations with circulating protein levels, with noncoding variants equally likely to increase or decrease protein levels (unlike coding variants which predominantly decrease function) [70]. Aggregate testing of noncoding variants identified 357 conditionally independent associated regions, 21% of which were not detectable by single-variant testing alone, highlighting the complementary value of burden approaches for noncoding variant interpretation.

The integration of UK Biobank and All of Us resources has created an unprecedented opportunity to advance gene-based burden testing for rare variant discovery. Standardized protocols for variant annotation, burden score construction, and statistical analysis enable robust identification of genes influencing complex traits and diseases. The cross-biobank validation framework enhances confidence in novel associations. Future directions will likely focus on improved functional interpretation of noncoding variants, development of more sophisticated burden testing methods that integrate common and rare variant signals, and expanded diverse representation to ensure equitable translation of genomic discoveries. These biobanks continue to serve as foundational resources for elucidating the genetic architecture of human disease and identifying novel therapeutic targets.

The journey from gene discovery to drug target validation represents a critical pathway in modern therapeutic development, particularly for rare genetic diseases. Gene-based burden testing for rare variants has emerged as a powerful statistical approach for identifying novel disease-associated genes by aggregating the effects of multiple rare genetic variants within a gene or genomic region [66]. This methodology addresses the fundamental challenge of limited statistical power when studying individual rare variants in small cohorts, enabling researchers to detect significant genetic associations that would otherwise remain obscured [66] [74]. The integration of these discoveries into validated drug targets requires a systematic framework that combines genomic evidence with functional validation and assessment of therapeutic potential, forming a continuum from initial genetic observation to target qualification for drug development programs [75] [76].

Application Notes: Integrating Burden Testing into the Therapeutic Development Pipeline

Key Statistical Methods for Gene-Based Burden Testing

Gene-based burden tests address the limitations of single-variant association analyses by collapsing multiple rare variants into a single test unit. Table 1 summarizes the primary statistical approaches utilized in rare variant association testing, each with distinct advantages and applications in gene discovery workflows [66].

Table 1: Statistical Methods for Rare Variant Burden Testing

Method	Type	Key Features	Advantages	Limitations
Burden Tests	Collapsing	Aggregates rare variants into a single score	High power when all variants affect trait in same direction	Loses power with bidirectional or non-causal variants
SKAT	Kernel-based	Models distribution of variant effects; allows covariate adjustment	Handles bidirectional effects; flexible modeling	Lower power if all variants affect in same direction
SKAT-O	Omnibus (Hybrid)	Combines SKAT and burden tests using Fisher's method	Balances power across different genetic architectures	May lose power when few trait-associated variants exist
C-alpha Test	Distribution-based	Tests for variability in effect direction among variants	Detects both risk-increasing and protective variants	Lower power when effects are unidirectional
ACAT	p-value Combination	Combines p-values using Cauchy distribution	Good power when few strong-effect variants are present	May underperform with many weak signals

From Genetic Association to Therapeutic Hypothesis: A Structured Workflow

The transition from a statistically significant gene-based burden test result to a therapeutic hypothesis requires multiple validation steps. First, genes identified through burden testing must be prioritized based on genetic evidence strength, with PathVar-identified single nucleotide variants (SNVs) demonstrating particular value in distinguishing case cohorts from controls, as evidenced in hemiplegic migraine research where this approach identified 20 significantly associated genes [66]. Second, biological plausibility must be established through pathway analysis and literature review, with novel gene discoveries such as SLC38A10, GCOM1, and NXPH2 in hemiplegic migraine requiring functional characterization to understand their potential roles in disease mechanisms [66]. Finally, practical considerations for therapeutic development must be assessed, including gene expression patterns across tissues, the nature of variant consequences (loss-of-function vs. gain-of-function), and potential safety concerns based on known gene functions in model organisms or human populations [75] [76].

Incorporating 3D Genomics for Enhanced Target Prioritization

Advanced genomic technologies now enable more confident linking of non-coding variants identified in burden studies to their target genes. Three-dimensional multi-omics approaches map the physical folding of the genome to identify long-range interactions between regulatory elements and genes, providing critical functional context for disease-associated variants [77]. This methodology addresses a fundamental limitation of traditional genomics, which often incorrectly assumes that variants affect the nearest gene in the linear DNA sequence - an approach that fails approximately 50% of the time according to Dr. Dan Turner of Enhanced Genomics [77]. By integrating genome folding data with chromatin accessibility and gene expression information, researchers can build a comprehensive regulatory network underlying disease pathogenesis, significantly enhancing confidence in target selection for therapeutic development.

Experimental Protocols

Protocol 1: Gene-Based Burden Testing Using Whole Exome Sequencing Data

Purpose and Scope

This protocol details a standardized workflow for conducting gene-based burden tests using whole exome sequencing (WES) data to identify genes enriched for rare variants associated with disease phenotypes. The methodology is adapted from recent research on hemiplegic migraine that successfully identified novel gene associations through this approach [66].

Materials and Equipment

• Whole exome sequencing data from cases and controls (minimum recommended sample size: 184 individuals) [66]
• Bioinformatics pipeline for variant calling (e.g., GATK best practices workflow)
• Variant annotation tools (e.g., PathVar, ANNOVAR, VEP)
• Population frequency databases (gnomAD, 1000 Genomes)
• Statistical software environment (R, Python) with burden testing packages (SKAT, SKAT-O, burden)
• High-performance computing resources for large-scale genetic analyses

Step-by-Step Procedure

• Variant Quality Control and Filtering

  • Process raw sequencing data through standard quality control metrics
  • Retain variants with call rate >95% and Hardy-Weinberg equilibrium p-value >1×10⁻⁶ in controls
  • Filter for rare variants with minor allele frequency (MAF) <0.01 in population databases
  • Apply additional quality filters based on sequencing depth, genotype quality, and mapping quality

• Variant Annotation and Prioritization

  • Annotate variants using PathVar or similar pipeline to identify pathogenic variants
  • Categorize variants by functional impact (missense, loss-of-function, splice-site)
  • Apply in silico prediction scores (SIFT, PolyPhen-2, CADD) to prioritize deleterious variants

• Gene-Based Collapsing Analysis

  • Aggregate qualifying rare variants within each gene for cases and controls separately
  • Implement multiple burden testing strategies: Strategy 1: PathVar-identified SNVs only Strategy 2: PathVar SNVs combined with missense and rare variants [66]
  • Adjust for population stratification using principal component analysis (PCA)

• Statistical Analysis and Significance Testing

  • Apply burden tests, SKAT, and SKAT-O to collapsed gene scores
  • Correct for multiple testing using Bonferroni correction or false discovery rate (FDR)
  • Define significance threshold at FDR <0.05
  • Validate associations in independent cohorts when available

Anticipated Results and Interpretation

The sequential application of multiple burden testing strategies typically yields a refined list of significantly associated genes. In the hemiplegic migraine study, Strategy 1 identified 20 significant genes, while the more stringent Strategy 2 refined this to 11 genes, demonstrating how integrative approaches reduce false positives while highlighting high-confidence associations [66]. Genes with burden test p-values below the significance threshold after multiple testing correction represent high-priority candidates for functional validation and therapeutic exploration.

Protocol 2: Functional Validation of Candidate Genes Using TAG Degradation

Purpose and Scope

This protocol describes a target validation methodology using TAG degradation technology to rapidly assess the phenotypic consequences of candidate gene knockdown in relevant cellular models. This approach provides critical functional evidence for genes identified through burden testing, establishing a direct link between genetic association and disease biology [78].

Materials and Equipment

• TAG Degradation Platform (dTAG, aTAG, or BromoTag systems) [78]
• Cell line relevant to disease pathology (primary cells or established lines)
• CRISPR-Cas9 system for locus-specific TAG knock-in
• TAG Degrader molecules (PROTACs, SNIPERs, PHOTACs, or uSMITE molecules)
• Simple Western system or conventional western blot equipment [78]
• Validated antibodies for target protein detection
• Cell culture equipment and appropriate media
• Phenotypic assay reagents (cell viability, functional assays)

Step-by-Step Procedure

• Generation of TAG Knock-In Cell Lines

  • Design CRISPR guide RNAs to fuse TAG domain to candidate gene via homologous recombination
  • Transfect cells with CRISPR components and TAG donor template
  • Select successfully edited clones using antibiotic resistance or fluorescence markers
  • Validate TAG fusion expression by western blot and Sanger sequencing

• TAG Degrader Treatment Optimization

  • Culture TAG knock-in cells in appropriate medium
  • Treat cells with titrated concentrations of TAG Degrader molecule (typical range: 1 nM-1 μM)
  • Include DMSO-only treated cells as negative controls
  • Harvest cells at multiple time points (e.g., 2, 4, 8, 12, 24 hours) to assess degradation kinetics

• Assessment of Target Knockdown

  • Lyse cells and quantify protein concentration
  • Analyze target protein levels using Simple Western or conventional western blot [78]
  • Determine DCâ‚…â‚€ (concentration achieving 50% degradation) and Dmax (maximum degradation)
  • Validate degradation specificity using non-TAGged control cells

• Phenotypic Characterization

  • Implement disease-relevant functional assays following target degradation
  • Assess cellular phenotypes at multiple time points post-degrader treatment
  • Include rescue experiments with degradation-resistant target variants when possible
  • Document phenotypic changes with appropriate statistical analysis

Anticipated Results and Interpretation

Successful target validation demonstrates a direct relationship between candidate gene degradation and relevant phenotypic changes. Dose-dependent reduction in target protein levels with corresponding DCâ‚…â‚€ values typically in the nanomolar range provides evidence of efficient degradation [78]. Phenotypic assays should reveal consequences consistent with the disease pathology and genetic association - for example, neuronal hyperexcitability for migraine-associated genes or ion transport defects for channelopathy genes. These functional data substantially increase confidence in the therapeutic relevance of candidate genes identified through burden testing.

Research Reagent Solutions

Table 2: Essential Research Tools for Gene Discovery and Target Validation

Category	Specific Product/Platform	Key Function	Application Notes
Bioinformatics Analysis	Geneyx Analysis Platform	NGS data analysis and interpretation	Accelerates diagnostic process from raw genetic data to clinical report; supports WGS, WES, gene panels [79]
Variant Prioritization	PathVar Pipeline	Identifies pathogenic variants from WES data	Critical for burden testing strategies; distinguishes case cohorts through distinctive SNV profiles [66]
Target Validation	TAG Degradation Platform (dTAG, aTAG, BromoTag)	Induces rapid degradation of target proteins	Enables functional validation of candidate genes; useful in cell culture and in vivo [78]
Protein Detection	Simple Western System	Automated protein separation and detection	Quantifies target protein knockdown (DCâ‚…â‚€, Dmax); integrates into degradation workflows [78]
3D Genomics	Enhanced Genomics Platform	Maps genome folding and regulatory interactions	Links non-coding variants to target genes; identifies causal regulatory networks [77]
Statistical Analysis	SKAT/SKAT-O R/Python Packages	Implements rare variant association tests	Provides flexible burden testing frameworks with covariate adjustment capabilities [66]

Workflow and Pathway Diagrams

From Gene Discovery to Target Validation Workflow

Statistical Framework for Burden Testing

TAG Degradation Mechanism

Independent Replication and Functional Validation of Burden Test Findings

Gene-based burden tests are a fundamental analytical method in rare-variant association studies, designed to elucidate the genetic architecture of complex diseases and traits. These tests operate by aggregating multiple rare genetic variants within a specific gene or genomic region into a single composite "burden" score, which is then tested for association with a phenotype of interest [26]. This approach is particularly powerful for detecting the collective influence of rare coding variations (typically with minor allele frequency, MAF < 1%) that individually have low statistical power due to their rarity but may exert significant effects on protein function and disease risk when considered jointly [66]. The core premise of burden testing is that multiple rare, deleterious variants within a functional unit (e.g., a gene) are more likely to be enriched in cases than controls and act in the same direction of effect on the trait [17].

Against the backdrop of large-scale exome and genome sequencing efforts in biobanks like the UK Biobank and All of Us, burden tests have successfully identified novel gene-trait associations for a wide spectrum of conditions, including metabolic diseases [80], migraine disorders [66], and quantitative blood biomarkers [17]. However, the initial statistical association represents only the first step in the gene discovery pipeline. Independent replication in separate cohorts and functional validation using experimental models are critical, sequential steps required to transform a statistical observation into a biologically and clinically meaningful finding [81]. This document outlines detailed protocols and application notes for confirming and validating burden test discoveries, providing a framework for researchers and drug development professionals to robustly characterize novel gene-disease relationships.

Statistical Replication Frameworks

Meta-Analysis of Rare-Variant Associations

Independent replication through meta-analysis is essential for verifying initial burden test findings. Meta-analysis combines summary statistics from multiple studies, dramatically increasing sample size and power to detect genuine rare-variant associations that may be underpowered in individual studies [82].

Table 1: Software Packages for Meta-Analysis of Rare-Variant Associations

Software Package	Statistical Methodology	Key Features	Input Requirements
RAREMETAL	Score-based burden and variance component tests	Efficient p-value computation, work with related samples	Single-variant score statistics, covariance matrices
MetaSKAT	Optimal sequence kernel association test (SKAT-O)	Combines burden and SKAT tests, handles bidirectional effects	Single-variant score statistics, covariance matrices
MASS	Burden and variance component tests	Flexible weighting schemes, user-friendly interface	Single-variant summary statistics
seqMeta	Burden tests and SKAT	R package, integrates with SeqArray data format	Single-variant summary statistics

The meta-analysis workflow involves several critical steps. First, participating studies calculate single-variant summary statistics, including effect sizes (beta coefficients), standard errors, and covariance matrices between variants [82]. For quantitative traits, appropriate transformations (e.g., inverse normal transformation) should be applied to ensure normality assumptions are met. Consortiums then combine these summary statistics using fixed-effects or random-effects models, with careful consideration of population structure and relatedness within studies. Finally, conditional analyses should be performed to remove "shadow" association signals driven by linkage disequilibrium with common variants, ensuring the rare-variant signal is genuinely independent [82].

The Sparse Burden Association Test (SBAT) for Replication

The Sparse Burden Association Test (SBAT) is a recently developed method that addresses key challenges in replicating burden test findings. Traditional burden testing approaches often calculate numerous correlated burden scores across different allele frequency thresholds and annotation classes, complicating multiple testing correction and interpretation [26]. SBAT addresses this by jointly testing a set of burden scores under the assumption that causal burden scores act in the same effect direction, using a non-negative least squares (NNLS) approach with positivity constraints on regression parameters [26].

The SBAT methodology can be implemented using the REGENIE software (v.3.0+) with the following protocol:

• Input Preparation: Format genotype data into a matrix of burden scores (X), where each column represents a distinct burden score derived from variants within specific AAF bins and annotation classes.
• Model Fitting: Fit the NNLS model: Y = Xβ + ε, where β ≥ 0, using the active set method to obtain estimates β̂ [26].
• Hypothesis Testing: Test H0: β = 0 versus H1: β ≥ 0 using the quadratic form test statistic T = β̂′V⁻¹β̂, where V = σ²(X′X)⁻¹ is the covariance matrix.
• p-value Calculation: Compute p-values from the null distribution, which is a mixture of chi-square distributions with weights dependent on the covariance between burden scores.
• Directionality Adjustment: Apply SBAT twice to both Y and -Y, then combine the two p-values using the Cauchy combination method to account for unknown effect direction [26].

Table 2: Comparison of Burden Test Methods for Replication Studies

Method	Key Principle	Advantages for Replication	Limitations
Standard Burden Test	Collapses variants into single score	Simple interpretation, high power when all variants are causal and same direction	Loses power with bidirectional effects or non-causal variants
SKAT	Variance-component test	Powerful for bidirectional variant effects, flexible modeling	Lower power when all variants affect trait in same direction
SKAT-O	Hybrid of burden and SKAT	Balances power across genetic architectures	May lose power with few trait-associated variants
SBAT	Joint testing with non-negativity constraints	Model selection sparsity, handles correlated burden scores	Requires effect direction consistency, newer method

Figure 1: Statistical replication workflow for burden test findings, integrating meta-analysis and sparse burden association testing.

Functional Validation Protocols

CRISPR-Based Functional Validation

For variants of unknown significance (VUS) identified through burden testing, CRISPR-Cas9 gene editing provides a direct method for functional validation in cellular models. This approach was successfully demonstrated in a Kleefstra syndrome case study involving the EHMT1 gene [83].

Protocol: CRISPR Validation of Gene Haploinsufficiency

• Guide RNA Design: Design and synthesize guide RNAs (gRNAs) targeting the exon or functional domain affected by the identified rare variants. For haploinsufficiency disorders, aim to introduce frameshift mutations that mimic predicted loss-of-function variants.
• Cell Line Selection: Select an appropriate cell model. HEK293T cells were used for Kleefstra syndrome validation [83], but disease-relevant cell types (e.g., neurons, iPSCs) are preferred when available.
• Transfection and Clonal Selection: Transfect cells with CRISPR-Cas9 components and perform high-throughput clonal selection to isolate isogenic cell lines with the desired edit. Validate edits using Sanger sequencing or next-generation sequencing.
• Transcriptomic Profiling: Extract total RNA from edited and control cells and perform RNA sequencing (RNA-seq) to identify differentially expressed genes (DEGs) and pathway alterations.
• Phenotypic Correlation: Compare the transcriptomic signature to the known disease phenotype. In the Kleefstra syndrome example, researchers identified changes in neural gene expression and chromosome 19/X suppression consistent with the clinical phenotype [83].

Multi-Omics Functional Assessment

Integrating multiple functional genomics approaches provides complementary evidence for variant pathogenicity, particularly when CRISPR validation is not feasible.

Protocol: Multi-Tiered Functional Genomics

• Transcriptomic Analysis: Perform RNA-seq on patient-derived cells (e.g., fibroblasts) to detect aberrant splicing, allele-specific expression, or gene expression changes. This approach increased diagnostic yield by 10% in mitochondrial disorders compared to exome sequencing alone [81].
• Biomarker Validation: For metabolic traits, measure relevant biomarkers in carrier plasma or cell culture media. For example, burden tests identifying ALPL associations with alkaline phosphatase should be followed by enzymatic activity assays [17].
• Protein Functional Assays: For missense variants, express mutant proteins and assess functional properties. For ion channel genes implicated in hemiplegic migraine, electrophysiology studies would be appropriate to characterize channel dysfunction [66].
• Pathway Analysis: Integrate omics data to identify disrupted biological pathways. Enrichment analysis of DEGs can reveal whether the gene's function aligns with the disease mechanism.

Figure 2: Functional validation pathways for burden test findings, combining gene editing and multi-omics approaches.

Application Notes: Case Studies in Complex Diseases

Cardiometabolic Traits

Large-scale whole-genome sequencing (WGS) studies in UK Biobank and All of Us have demonstrated the utility of burden testing for identifying novel genes for cardiometabolic diseases. Recent WGS analyses of ~490,000 individuals identified protein-truncating variants (PTVs) in IRS2 with substantial effects on type 2 diabetes risk (OR = 6.4) and chronic kidney disease independent of diabetes status [80]. Similarly, PTVs in RIF1 were associated with body mass index (effect = 2.66 kg/m²) [80]. These findings highlight several important principles for replication and validation:

• Sensitivity Analyses: Adjust for polygenic risk scores and regional common variants to ensure rare-variant signals are independent. For IRS2, associations persisted after these adjustments, supporting a genuine allelic series [80].
• Leave-One-Out Analysis: Verify that gene-level associations are not driven by a single rare variant by iteratively removing each variant and recalculating associations.
• Cross-Ancestry Validation: Test associations in diverse populations. WGS provided stronger associations than WES, partly due to increased ancestral diversity in analyzed samples [80].

Neurological Disorders

In hemiplegic migraine (HM), a rare neurovascular disorder, burden testing of whole-exome sequencing data from 184 patients identified novel genes including SLC38A10, GCOM1, and NXPH2, expanding the genetic architecture beyond known ion channel genes [66]. The validation strategy included:

• Variant Prioritization: Using PathVar, a bioinformatics pipeline designed to identify pathogenic variants through functional annotation.
• Burden Test Implementation: Applying multiple association strategies: (1) PathVar-identified SNVs only and (2) PathVar SNVs combined with missense and rare variants.
• Biological Plausibility Assessment: Evaluating whether newly associated genes fit the known disease mechanism of neuronal excitability and cortical spreading depression.

Blood Biomarkers

For quantitative blood biomarkers, gene-based burden scores have identified both positive and negative effect directions. For example, rare damaging variants in ALPL decrease alkaline phosphatase levels, while variants in LDLR increase LDL direct measurement [17]. Validation approaches included:

• Directionality Consistency: Confirming that effect directions align with biological expectations (e.g., damaging variants in protein-coding genes decrease corresponding biomarker levels).
• Multi-Method Comparison: Comparing results across different association methods (SAIGE-GENE, SKAT-O, Fisher's exact test) to ensure robustness.
• Clinical Correlation: Relating genetic findings to known disease mechanisms (e.g., LDLR mutations in familial hypercholesterolemia) [17].

The Scientist's Toolkit

Table 3: Essential Research Reagents and Platforms for Burden Test Validation

Category	Specific Tools/Reagents	Function in Validation	Application Examples
Statistical Genetics Software	REGENIE (SBAT implementation) [26]	Joint testing of correlated burden scores with sparsity	Replication of multi-threshold burden tests
	RAREMETAL, MetaSKAT [82]	Meta-analysis of rare-variant associations	Combining summary statistics across cohorts
Gene Editing Tools	CRISPR-Cas9 with gRNA constructs [83]	Introduce specific variants in cellular models	Functional testing of VUS in disease genes
Omics Technologies	RNA-sequencing platforms	Transcriptome profiling of edited cells/carriers	Identify pathway disruptions (e.g., neural genes in Kleefstra syndrome [83])
	REVEL (Rare Exome Variant Ensemble Learner) [80]	In silico prediction of missense variant pathogenicity	Variant prioritization for functional testing
Biomarker Assays	ELISA/MS-based protein quantification	Measure biomarker levels in carrier samples	Validate ALPL-alkaline phosphatase association [17]
Cell Models	HEK293T, patient-derived fibroblasts, iPSCs	Cellular context for functional studies	Functional validation for metabolic/neurological disorders [81] [83]

Conclusion

Gene-based burden testing has matured into an indispensable tool for deciphering the role of rare variants in human disease, successfully bridging a critical gap in complex trait genetics. The development of robust, scalable methods that handle technical artifacts and case-control imbalance, coupled with efficient meta-analysis frameworks, has solidified the reliability of this approach. Looking forward, the integration of burden test results with other data types—such as common variant polygenic risk scores and real-world evidence—holds immense promise for creating more comprehensive genetic risk models. For drug development, these tests are already providing crucial human genetic evidence to prioritize and de-risk therapeutic targets, as exemplified by frameworks predicting on-target side effects. The continued expansion of biobank resources and the advent of advanced, statistically powerful tests will further accelerate the translation of genetic discoveries into clinical applications and novel treatments.

References

• 1. Rare-Variant Association Analysis: Study Designs and ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4085641/]
• 2. Unique roles of rare variants in the genetics of complex ... [https://www.nature.com/articles/s10038-020-00845-2]
• 3. Rare disease gene association discovery in the ... [https://www.nature.com/articles/s41586-025-08623-w]
• 4. Gene-Based Burden Testing of Rare Variants in ... [https://pubmed.ncbi.nlm.nih.gov/40725463/]
• 5. COBT: A gene-based rare variant burden test for case-only ... [https://www.medrxiv.org/content/10.1101/2025.07.16.25331553v1]
• 6. Prioritizing disease-related rare variants by integrating gene ... [https://pubmed.ncbi.nlm.nih.gov/39348415/]
• 7. Gene-Based Burden Testing of Rare Variants in ... [https://www.mdpi.com/2073-4425/16/7/807]
• 8. A comprehensive evaluation of collapsing methods using ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC4164031/]
• 9. A new approach for rare variation collapsing on functional ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6499321/]
• 10. Gene-based burden tests of rare germline variants identify ... [https://pubmed.ncbi.nlm.nih.gov/39472694/]
• 11. Genome-wide rare variant analysis for thousands of ... [https://www.nature.com/articles/s41467-020-14288-y]
• 12. Article A case–control collapsing analysis identifies retinal ... [https://www.sciencedirect.com/science/article/pii/S1098360021044981]
• 13. Strategies for aggregating gene expression data [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-12-322]
• 14. Burden Testing of Rare Variants Identified through Exome ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6174288/]
• 15. An optimized variant prioritization process for rare disease ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12539062/]
• 16. Baylor Genetics Presents New Data on Clinical and ... [https://www.baylorgenetics.com/news/baylor-genetics-presents-new-data-on-clinical-and-diagnostic-utility-of-rna-sequencing-for-rare-disease-at-ashg-2025-annual-meeting/]
• 17. Gene-based burden scores identify rare variant associations ... [https://bmcgenomdata.biomedcentral.com/articles/10.1186/s12863-023-01155-0]
• 18. Scalable and accurate rare variant meta-analysis with Meta ... [https://www.nature.com/articles/s41588-025-02403-y]
• 19. Statistical Analysis Strategies for Association Studies ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3743540/]
• 20. Recent advances and challenges of rare variant ... [https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2022.1014947/full]
• 21. Rare-variant association studies: When are aggregation ... [https://www.sciencedirect.com/science/article/pii/S0002929725002721]
• 22. Winner's curse in rare variant analysis: effect size ... [https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2025.1416673/full]
• 23. Real world scenarios in rare variant association analysis [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2591-6]
• 24. Rare variant association studies: Significance, methods ... [https://www.sciencedirect.com/science/article/pii/S106345842401505X]
• 25. Optimal tests for rare variant effects in sequencing association ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC3440237/]
• 26. Joint testing of rare variant burden scores using non- ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11480795/]
• 27. SAIGE-GENE+ improves the efficiency and accuracy of set- ... [https://www.nature.com/articles/s41588-022-01178-w]
• 28. Genome-wide prediction of disease variant effects with a ... [https://www.nature.com/articles/s41588-023-01465-0]
• 29. Fine-tuning protein language models to understand the ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12166733/]
• 30. Article Landscapes of missense variant impact for human ... [https://www.sciencedirect.com/science/article/pii/S0002929725003519]
• 31. Evaluation of GENESIS, SAIGE , REGENIE and fastGWA-GLMM for... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9544087/]
• 32. SAIGE-GENE+ improves the efficiency and accuracy of set- ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9534766/]
• 33. Evaluation of GENESIS, SAIGE, REGENIE and fastGWA ... [https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2022.897210/full]
• 34. - SAIGE + improves the efficiency and accuracy... | Nature Genetics GENE [https://www.nature.com/articles/s41588-022-01178-w?error=cookies_not_supported&code=5ed50ffb-165e-4c52-9574-3976e8c64676]
• 35. Burden module for mrcepid-runassociationtesting [https://github.com/mrcepid-rap/burden]
• 36. Population-scale gene-based analysis of whole-genome ... [https://www.nature.com/articles/s41588-025-02364-2]
• 37. Identifying genetic risk factors for serious adverse drug ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC2763923/]
• 38. A gene based combination test using GWAS summary data [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-022-05114-x]
• 39. Genetic ancestry in population pharmacogenomics ... [https://www.sciencedirect.com/science/article/pii/S2589004224021412]
• 40. Controlling for human population stratification in rare variant ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8463695/]
• 41. Computationally efficient meta-analysis of gene-based ... [https://www.nature.com/articles/s41588-025-02390-0]
• 42. Composite quantile regression approach to batch effect ... [https://www.frontiersin.org/journals/microbiology/articles/10.3389/fmicb.2025.1484183/full]
• 43. A Fast and Accurate Algorithm to Test for Binary Phenotypes ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5501775/]
• 44. Efficiently controlling for case-control imbalance and ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC6119127/]
• 45. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies - PubMed [https://pubmed.ncbi.nlm.nih.gov/30104761/]
• 46. COBT: A gene-based rare variant burden test for case-only ... [https://www.medrxiv.org/content/10.1101/2025.07.16.25331553v1.full-text]
• 47. Leveraging functional annotations to map rare variants ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11643288/]
• 48. Leveraging functional annotations to map rare variants ... [https://www.sciencedirect.com/science/article/pii/S0002929725002861]
• 49. Effective filtering strategies to improve data quality from ... [https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-15-125]
• 50. Gene-based burden tests of rare germline variants identify ... [https://www.nature.com/articles/s41588-024-01966-6]
• 51. Integrating variant functional annotation scores have varied ... [https://www.nature.com/articles/s41598-022-14924-1]
• 52. Scalable and Accurate Meta-Analysis for Rare Variants [https://www.medrxiv.org/content/10.1101/2024.09.17.24313855v1]
• 53. Joint testing of rare variant burden scores using non- ... [https://pubmed.ncbi.nlm.nih.gov/39366334/]
• 54. Article Joint testing of rare variant burden scores using non ... [https://www.sciencedirect.com/science/article/pii/S0002929724003070]
• 55. A Comparison of Methods for Gene-Based Testing That ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC9117705/]
• 56. an evaluation of power and type I error of gene-based tests ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC5133510/]
• 57. Gene-based association tests using GWAS summary ... [https://www.nature.com/articles/s41598-022-07465-0]
• 58. Comparison of methods for multivariate gene-based ... [https://www.nature.com/articles/s41431-018-0327-8]
• 59. A Rapid and Powerful Gene-Based Association Test Using ... [https://www.sciencedirect.com/science/article/pii/S0002929711000498]
• 60. Specificity, length and luck drive gene rankings in ... [https://www.nature.com/articles/s41586-025-09703-7]
• 61. Combining rare and common genetic variants improves ... [https://www.sciencedirect.com/science/article/pii/S2949774424009725]
• 62. Contributions of Common, Rare, and Somatic Genetic ... [https://pubmed.ncbi.nlm.nih.gov/41060636/]
• 63. New AI model could speed rare disease diagnosis [https://www.broadinstitute.org/news/new-ai-model-could-speed-rare-disease-diagnosis]
• 64. New PopEVE AI Model Speeds Up Rare Disease Discovery [https://oncodaily.com/industry/popeve]
• 65. Rare variants: data types and analysis strategies - PMC [https://pmc.ncbi.nlm.nih.gov/articles/PMC8267277/]
• 66. Gene-Based Burden Testing of Rare Variants in Hemiplegic ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12294318/]
• 67. Lessons from national biobank projects utilizing whole ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC11887102/]
• 68. Weighted burden analysis of rare coding variants in ... [https://www.nature.com/articles/s10038-024-01235-8]
• 69. Rare coding variant analysis for human diseases across ... [https://pubmed.ncbi.nlm.nih.gov/39210047/]
• 70. Whole-genome sequencing analysis identifies rare, large ... [https://www.nature.com/articles/s41588-025-02095-4]
• 71. Joint testing of rare variant burden scores using non- ... [https://www.ukbiobank.ac.uk/publications/joint-testing-of-rare-variant-burden-scores-using-non-negative-least-squares/]
• 72. Systematic single-variant and gene-based association ... [https://www.sciencedirect.com/science/article/pii/S2666979X22001100]
• 73. All of Us Research Program year in review: 2024 [https://www.sciencedirect.com/science/article/pii/S0002929725003131]
• 74. Strategies to Uplift Novel Mendelian Gene Discovery for ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC8248347/]
• 75. Target Validation - NCBI - NIH [https://www.ncbi.nlm.nih.gov/books/NBK195039/]
• 76. Improving target assessment in biomedical research [https://www.nature.com/articles/s41573-020-0087-3]
• 77. Mapping the genome in 3D to reveal new drug targets [https://www.drugtargetreview.com/article/190646/mapping-the-genome-in-3d-to-reveal-new-drug-targets/]
• 78. Target Validation and Exploration [https://www.bio-techne.com/research-areas/targeted-protein-degradation/target-validation-exploration]
• 79. Geneyx: NGS Data Analysis | Clinical Genome Research [https://geneyx.com/]
• 80. Population-scale gene-based analysis of whole-genome ... [https://pmc.ncbi.nlm.nih.gov/articles/PMC12513836/]
• 81. functional validation of genetic variants - PMC - PubMed Central [https://pmc.ncbi.nlm.nih.gov/articles/PMC5959958/]
• 82. Meta-analysis for Discovering Rare-Variant Associations [https://www.sciencedirect.com/science/article/pii/S000292971500186X]
• 83. Functional validation of variants of unknown significance ... [https://www.sciencedirect.com/science/article/pii/S0378111922001068]