Computationally Efficient Whole-Genome Signal Region Detection for Quantitative and Binary Traits

Understanding the genetic underpinnings of human diseases and traits remains a central focus of genetic research. Genome-Wide Association Studies (GWAS) have been instrumental in exploring the genetic architecture of complex diseases and traits over the past decade (Visscher et al., 2017). By using array-based technologies, GWAS analyzes millions of single nucleotide polymorphisms (SNPs) across the genome to identify those associated with specific traits or disease outcomes. While GWAS has successfully identified thousands of common genetic variants linked to disease susceptibility, these common variants explain only a small fraction of heritability, which is often referred to as the “missing heritability problem”(Manolio et al., 2009). The majority of genetic variants in the human genome are rare, and these rare variants are believed to contribute significantly to the unexplained heritability. However, the classical GWAS approach which focuses on single-SNP-based analysis has very limited power for analyzing rare variants, as the effect of a single SNP can be too small to be detected. (Liu et al., 2010; Bakshi et al., 2016).

To address this limitation, whole genome sequencing (WGS) studies are being conducted to identify rare variants associated with disease susceptibility. Progress has been made in developing set-based association methods, which test multiple variants jointly by aggregating their effects within defined genomic regions. These methods include burden tests (Morgenthaler & Thilly, 2007; Madsen & Browning, 2009), the Sequence Kernel Association Test (SKAT) (Wu et al., 2011), and STAAR which incorporates variant functional annotations to enhance detection power (Li, Li, Zhou et al., 2022). The STAAR-O test extends the STAAR framework as an omnibus test by combining multiple annotation-weighted methods into a single unified test. A common challenge in these approaches is defining regions for variant sets, especially in non-coding or intergenic regions without clear functional boundaries. STAAR addresses this issue by applying gene-centric analysis to well-defined gene-associated regions and fixed sliding window-based analysis to regions without clear boundaries. Alternative approaches include the scan statistic (Naus, 1982), which systematically searches the human genome using a fixed window size. Mean-based scan statistic methods are later proposed for DNA copy number analysis, allowing the use of multiple window sizes in settings closely related to change-point detection problems (Jeng et al., 2010; Olshen et al., 2004; Zhang et al., 2010). Other notable frameworks, such as those based on the knockoff (i.e.KnockoffScreen) (He et al., 2021a, b), conducts genome-wide set-based analyses to identify signal regions while mitigating the effects of correlation confounding. However, fixed-window approaches can result in a loss of power because the sizes of signal regions can vary across the genome.

Further advancements have been made to enable signal region detection with dynamic window sizes. Li, Li, Zhou et al. (2020) introduced SCANG which combines the scan algorithm proposed by Jeng et al. (2010) with burden tests, SKAT, and the omnibus test to continuously scan the genome. However, this approach lacks comprehensive theoretical and empirical analysis of false discoveries and shows limited power when functional annotations are not incorporated (Li, Li, Zhou et al., 2022). Building on the scan algorithm, Li, Liu & Lin (2020) developed the quadratic scan statistic (QSCAN), which aggregates information across intervals of varying sizes. QSCAN provides theoretical guarantees for controlling false discoveries and has demonstrated strong detection performance, particularly when signal regions contain both causal and neutral variants. However, scan-based methods require calculating test statistics for many candidate intervals and applying a fixed threshold for selection, which is computationally intensive and tends to be conservative. SCANG-STAAR (Li, Li, Zhou et al., 2022) extends SCANG’s dynamic window scanning by incorporating multiple functional annotations via STAAR to boost the power of detecting rare variant associations. While functional annotations can enhance power and interpretability, their effectiveness relies heavily on the quality and relevance of the annotations, which may introduce biases, computational challenges, and increased resource demand.

Beyond the scan algorithm, Zhang, Wang & Yao (2023) proposed the binary and re-search (BiRS) algorithm for detecting signal regions with dynamic window sizes in whole genome sequencing (WGS) studies. The BiRS algorithm iteratively splits identified regions until a minimum size is reached, providing theoretical guarantees for FWER and FDR, which is shown more computationally efficient than scan-based algorithms and demonstrates superior power compared to QSCAN and KnockoffScreen. Impressively, BiRS surpasses SCANG-STAAR in detecting moderate or weak signals without requiring functional annotations. The combination of computational speed, increased detection power and theoretical rigor make BiRS as a significant advancement in signal region detection for WGS studies. The BiRS algorithm was combined with the DCF two-sample mean test (Xue & Yao, 2020) for direct application to binary traits. However, it has not yet been extended to correct for multi-dimensional covariates or to test for non-binary outcomes.

In this paper, we generalize the BiRS algorithm to a distributed version (dBiRS) to combine with a new infinity-norm test statistic based on the summary statistics obtained from a generalized linear model. As the proposed test statistic fundamentally differs from the DCF two-sample mean test, this extension is non-trivial. We develope new theoretical guarantees to ensure the size control of dBiRS while preserving the detection accuracy of the original BiRS. The dBiRS algorithm operates in two main stages, combining results from local and central machines to support parallel computing and maintain computational efficiency, similar to distributed learning frameworks (Cai & Wei, 2022). In the first stage, BiRS is applied within genomic blocks, where local machines detect signal regions by calculating test statistics and thresholds within each block. Instead of simply aggregating block-wise results using a global threshold, the detected signal regions, along with their corresponding test statistics and thresholds, are transferred to a central machine. In the second stage, the central machine evaluates the significance of each block by performing a second layer of BiRS based on block-wise infinity-norm test statistics. Finally, a new threshold is constructed by multiplier bootstrap to reassess the signal regions within the significant blocks. We have shown theoretically that the dBiRS algorithm is able to consistently identifies true signal regions under more general alternative structures and in the presence of model misspecification under the alternative hypothesis. Simulated studies also demonstrate that dBiRS is more accurate and robust than the state-of-the-art KnockoffScreen and scan procedures.

Finally, we apply the BiRS algorithm to analyze whole exome sequence (WES) data from the UK Biobank, aiming to identify signal regions for intelligence and prospective memory. The dBiRS algorthm identified signal regions involving 327 genes, including 84 of which were previously reported to be assocaited with intelligence. Notable discoveries include both common and rare variants linked to cognitive functions, such as those in COL16A1, CRTC2, and PTPRF. Variants in genes like CRTC2, BRWD1, and TOP2B were also identified, highlighting their roles in endothelial function, immune system regulation, and neuronal development, respectively. In addition, the analysis revealed rare variants in 22 genes not previously associated with intelligence. Some of these genes are linked to Alzheimer’s disease, brain connectivity, and neuronal development, providing novel insights into the genetic basis of cognitive traits. For prospective memory, the study identified eight novel genes, including OMA1 and CNGB3, which are associated with neuroimaging measurements and cognitive functions.

The rest of the article is organized as follows. In Section 2, we introduce the testing procedure under GLM and describe the proposed dBiRS algorithm. We conduct comprehensive simulations in Section 3 to demonstrate that the proposed method enjoys preferably numerical performance compared with existing approaches. In Section 4, we apply the dBiRS algorithm to conduct Whole Exome Sequencing (WES) analyses on intelligence and prospective memory, aiming to deepen our understanding of cognitive aging and uncover genetic factors contributing to the risk of neurodegenerative disorders. We conclude the article with a discussion in Section 5, while the theoretical properties of the proposed algorithm, including size control and detection consistency, are deferred to the Appendix. Technical assumptions, lemmas, and proofs of the theoretical results, along with additional simulation and application results, are provided in the Supplementary Material.

We conduct simulation studies to compare the proposed dBiRS method with the state-of-art Q-SCAN and KnockoffScreen procedure. We generate sequence data of European ancestry from 10,000 chromosomes across 5-megabase (Mb) regions using the calibrated coalescent model (Sabeti & Schaffner, 2014), and the total number of variants is 349,640. We evaluate the family-wise error rate (FWER), false discovery rate (FDR), signal region detection rate (DR), and true positive rate (TPR) for both continuous and binary phenotypes.

For the evaluation of FWER, the continuous phenotypes are generated using the model:

where $X_{1}$ is a continuous covariate sampled from a standard normal distribution, $X_{1}\sim\mathcal{N}(0,1)$, and $X_{2}$ is a dichotomous covariate that takes values $0$ and $1$ with equal probability. The random noise $\varepsilon$ is generated from a standard normal distribution, $\varepsilon\sim\mathcal{N}(0,1)$. The dichotomous phenotypes are generated using the following logistic regression model:

where $\mathbb{P}(Y=1)=(\mathbb{P}(Y_{1}=1),\dots,\mathbb{P}(Y_{n}=1))^{\top}$. We perform dBiRS and Q-SCAN analyses based on 1,000 Monte Carlo runs under the linear and logistic models to compute FWERs. For both dBiRS and Q-SCAN, the number of bootstrap iterations is set to 1,000. We do not calculate the FWER for KnockoffScreen because it only controls the FDR (He et al., 2021a). The empirical FWERs estimated for dBiRS and Q-SCAN are presented in Table 1 for significance levels $\alpha=0.01$ and $\alpha=0.05$. The FWERs of Q-SCAN and dBiRS are accurate at both significance levels, demonstrating that both methods effectively control the FWER.

We next evaluate the detection accuracy of dBiRS and compare its performance with Q-SCAN and KnockoffScreen. We randomly select four 5kb causal windows, denoted as $I_{i}$, $i=1,\dots,4$. We generate continuous and dichotomous phenotypes by

and

where $\boldsymbol{G}_{I_{1}},\dots,\boldsymbol{G}_{I_{4}}$ represent the genotypes of the causal windows, and $\beta_{I_{1}},\dots,\beta_{I_{4}}$ are the corresponding effect sizes. In each causal window, $10\%$ of the variants are randomly designated as causal, and each causal variant is assigned an effect size as a decreasing function of the minor allele frequency (MAF), i.e., $\left\lvert\beta\right\rvert=c\left\lvert\log_{10}(\mathrm{MAF})\right\rvert$. The parameter $c$ is set to $c\in\left\{0.12,0.15\right\}$ for continuous outcomes and $c\in\left\{0.25,0.30\right\}$ for dichotomous outcomes. The sign of $\beta$ is randomly assigned with $50\%$ of the values being positive and $50\%$ negative.

We evaluate the DR,TPR and $h$ kilobase (kb) FDR for three methods based on 100 Monte Carlo runs across different values of $c$ (Li, Liu & Lin, 2022; He et al., 2021a). Specifically, we denote the true signal regions as $I_{1}^{*},\dots,I^{*}_{\ell}$ and the estimated signal regions as $\hat{I}_{1},\dots,\hat{I}_{\hat{\ell}}$. Let $I^{*}=\cup_{i=1}^{\ell}I^{*}_{i}$ represent the union of the true signal regions, and $\hat{I}=\cup_{i=1}^{\hat{\ell}}\hat{I}_{i}$ represent the union of the estimated signal regions. The DR and TPR are defined as

Following He et al. (2021a), we define the FDR($h$) for $h\in\{25,50,75\}$ as

to account for the spill-over effect of the LD structure, where $d(j,I^{*})=\min_{i\in I^{*}}\left\lvert i-j\right\rvert$ is the minimum distance between point $j$ and true signal regions. The DR represents the proportion of true signal regions that are detected and measures the ability to identify signal regions. In contrast, the TPR and FDR($h$) evaluate how well the true signal regions are recovered by assessing the similarity or distance between the true and identified signal regions. The DR, TPR, and FDR($h$) of the three methods for continuous and dichotomous phenotypes are presented in Table 2. Additionally, we provide the standard deviations (SDs) for DR, TPR, and FDR($h$) in Table 3.

Table 2 shows that dBiRS achieves the highest DRs and TPRs with the lowest FDRs across all signal strengths $c$ and for all values of the distance parameter $h=25,50,75$ in both continuous and dichotomous phenotypes. As the value of $h$ increases, the FDR of all methods decreases, and the FDR of dBiRS and Q-SCAN falls below 0.05 when $h=75$ kb. Q-SCAN outperforms KnockoffScreen in recovering true signal regions by achieving higher TPR and lower FDR. However, it identifies a smaller proportion of true signal regions (i.e., a lower DR) compared to KnockoffScreen. KnockoffScreen employs LD-pruning with a 0.75 correlation threshold to remove highly correlated variants before conducting the analysis. It identifies a larger number of signal regions, which increases the chance of overlap with true signal regions and thereby boosting the DR. However, the signal regions identified by KnockoffScreen have lower coverage of true signal regions due to both the LD-pruning process and a higher number of false discoveries. Additionally, Table 3 shows that dBiRS has the lowest standard deviations in all cases, indicating that dBiRS is the most stable method.

To further illustrate the recovery of true signal regions across different methods, we present the selection probabilities for all variants, where the selection probability of a variant $j$ is defined as the proportion of replications in which it is identified as a signal variant. We present the selection probabilities under the continuous phenotype and the dichotomous phenotype in Figure 1 and Figure 2, respectively.

While all procedures exhibit relatively high selection probabilities around the four true signal regions, dBiRS demonstrates the greatest accuracy and stability, achieving selection probabilities exceeding 0.95 in all true signal regions across all settings. When signals are relatively weak (i.e., $c=0.12$ for the continuous trait in Figure 1A and $c=0.25$ for the binary trait in Figure 2A), Q-SCAN shows relatively low selection probabilities in the leftmost true signal region (i.e., selection probability $<0.95$) and identifies some variants in non-signal regions at the beginning of the genome. This suggests that the Q-SCAN procedure is more prone to producing misleading results across replications, which is consistent with results that show a higher FDR compared to dBiRS. KnockoffScreen demonstrates the lowest selection probabilities in the four true signal regions and identifies only one variant with a selection probability greater than 0.95 (Figure 1A and Figure 2A). When the signal is stronger (i.e., $c=0.15$ for the continuous trait in Figure 1B and $c=0.30$ for the binary trait in Figure 2B), KnockoffScreen exhibits substantially higher selection probabilities in non-signal regions compared to both Q-SCAN and dBiRS (Figure 1B). The spread of selection probabilities across non-signal regions highlights a lack of stability for KnockoffScreen and indicates its tendency to identify false positives for spurious variants. Based on these results, we conclude that the dBiRS method outperforms the other procedures in terms of both detection accuracy and robustness across different settings.

Extensive research has identified common genetic variants associated with cognitive health, but the roles of rare genetic variants, which often have more subtle and complex effects, remain poorly understood. In this study, we apply the dBiRS algorithm to conduct Whole Exome Sequencing (WES) analyses on the core cognitive phenotypes: fluid intelligence (Field ID 20016) and prospective memory (Field ID 20018), using data from the UK Biobank (https://biobank.ctsu.ox.ac.uk/crystal/index.cgi). The WES data used in this study is derived from the final exome release in PLINK format (Field ID 23158) from the UK Biobank. Fluid intelligence reflects problem-solving and reasoning abilities, while prospective memory measures the capacity to remember and execute planned actions, a critical daily-life function that often declines with age. These phenotypes are closely linked to aging and neurodegeneration. By analyzing these traits, we aim to uncover rare variants that may confer either risk or resilience to cognitive functions, providing insights into the genetic and biological mechanisms underlying cognitive impairment. This approach has the potential to deepen our understanding of cognitive aging and uncover genetic factors contributing to the risk of neurodegenerative disorders such as Alzheimer’s disease (AD).

We implemented a standard quality control procedure to ensure the integrity of the dataset before analysis. First, we excluded individuals flagged as outliers by the UK Biobank based on genotyping missingness rates or heterogeneity, as well as those whose genotypically inferred sex did not align with their self-reported sex. To address population stratification, we utilized principal component analysis provided by the UK Biobank and excluded individuals identified as non-European. Specifically, we removed individuals whose values for either of the first two principal components deviated by more than five standard deviations from the mean. We also excluded participants who self-reported an ethnicity other than European. Furthermore, individuals with more than $5\%$ missing genotype data across variants passing UK Biobank’s quality control were removed. For variant-level quality control, we retained only biallelic autosomal variants assayed by both genotyping arrays used by the UK Biobank. Variants that failed UK Biobank quality control in any genotyping batch were excluded. Additionally, we removed variants with a Hardy-Weinberg equilibrium (HWE) $p$-value below $10^{-50}$ or with a minor allele count (MAC) $\leq 1$. This stringent filtering resulted in 13,681,006 variants across 22 chromosomes in the WES dataset, 154,785 samples for fluid intelligence analysis, and 155,448 samples for prospective memory analysis. After the quality control procedure, the dBiRS analysis was performed using summary statistics derived from a generalized linear model. The model adjusted for covariates, including sex, age, assessment center, and the top five genomic principal components provided by the UK Biobank.

Figure 1 illustrates the distribution of functional consequences of variants identified by dBiRS, with a substantial number of variants located in exonic, intronic, and UTR regions. Variants in the exonic region may directly alter protein structure and function. Variants in the 3’ UTR (UTR3) are likely involved in post-transcriptional regulation, including mRNA stability and microRNA binding. Intronic variants could disrupt splicing mechanisms or long-range regulatory elements, further impacting gene function. While exonic variants are commonly the focus of WES studies, our findings also underscore the important role of rare variants in non-coding and regulatory regions, such as UTRs and intronic regions, which have historically been understudied in the context of intelligence and prospective memory.

For the study on intelligence, dBiRS identified signal regions encompassing 327 genes, including 84 genome-wide significant genes previously reported for their association with intelligence. Additionally, there are 11 genome-wide significant genes overlapping with those associated with general cognitive function, 30 genes overlapping with those implicated in schizophrenia, 65 genes associated with educational attainment, and 35 genes overlapping with autism spectrum disorder. Based on the GWAS Catalog, 33 genome-wide significant common variants were also identified in our WES study. Notable examples include SNP rs2271928 in COL16A1, SNP rs11264680 in CRTC2, and SNP rs539096 in PTPRF, among others.

The dBiRS method identifies numerous additional rare variant associations in the exonic, intronic, and 3’ UTR regions of genes that are previously reported to be associated with intelligence or cognitive performance. For example, a previous study highlights the significance of CRTC2 in endothelial function, which is essential for maintaining vascular physiology in the brain (Kanki et al., 2020). While prior genome-wide studies have identified common variants in CRTC2 associated with cognitive traits (Savage et al., 2018), the dBiRS method uncovers 4 rare variants in the exonic region of CRTC2, which may directly alter protein function, 3 rare variants in the intronic region that could influence gene regulation, and 4 rare variants in the 3’ UTR region, potentially affecting mRNA stability and translation efficiency. Another example is BRWD1 that is previously identified as a significant gene associated with intelligence (Savage et al., 2018). In our study, we identified rare variants in its exonic, intronic, and 3’ UTR regions. Research has shown that BRWD1 plays a role in establishing epigenetic states during B cell development, which is essential for proper immune function (Mandal et al., 2018). Further studies are needed to clarify its specific contributions to neural development and cognition.

The dBiRS method identified rare variants in 22 novel genes that have not been previously reported in the GWAS Catalog for intelligence. Among these 22 genes, 14 were significant in the gene-based test for general cognitive ability (Davies et al., 2018). Of the remaining 9 genes, SZT2 has been linked to common variants associated with Alzheimer’s disease and specific brain regions (Herold et al., 2016); ITIH4 has been associated with schizophrenia (Ripke et al., 2011); and NDUFA6 has been linked to brain connectivity measurements (Zhang, Pan, Huang, Liang, Li, Zhang & Zhang, 2023). Genes SLC39A8, TOP2B, NKIRAS1, GLTBD1, and TPM3 are new findings that have not been previously reported in relation to cognitive performance measurements. Among these genes, TOP2B plays a critical role in neuronal development by regulating gene expression during brain development (Tiwari et al., 2012). Proper functioning of TOP2B is essential for neuronal connectivity and synaptic plasticity, both of which underlie learning and memory (Madabhushi et al., 2015). Disruptions in TOP2B have been linked to neurodevelopmental disorders, which could implicate it in intelligence (King et al., 2013).

For the study on prospective memory, we identified eight novel genes associated with prospective memory, a key component of cognitive function. Among these, OMA1 stands out due to its established links with multiple neuroimaging measurements, including white matter lesion progression and brain column structure, underscoring its potential role in cognitive processes (Wang et al., 2020). Notably, our analysis uncovered 10 rare variants in the 3’ UTR, 6 rare variants in the exonic region, 7 rare variants in the intronic region, and 1 common variant in the 5’ UTR of OMA1. These findings suggest that OMA1 may influence prospective memory through diverse mechanisms, such as protein function modulation, gene regulation, and mRNA stability. Similarly, we identified variants in the intronic, exonic, and 3’ UTR regions of CNGB3, further emphasizing the importance of regulatory regions in prospective memory. CNGB3 is known to be associated with educational attainment (Okbay et al., 2022), cognitive function (Lee et al., 2018) and anxiety (Meier et al., 2019), which supports its potential role in cognitive traits. Additionally, we discovered HIGD1B, previously reported to be associated with educational attainment (Pasman et al., 2022), and SFMBT2 which has been linked to anxiety and stress-related disorders (Hollins & Cairns, 2016). Other novel findings include the genes NT5C3B, NOSTRIN, and APOD, which represent unexplored candidates for prospective memory and related cognitive traits. These discoveries underscore the value of our approach in uncovering novel genetic mechanisms underlying complex cognitive traits and highlight new directions for future research into the genetic architecture of prospective memory.

We developed a computationally efficient distributed Binary and Re-Search (dBiRS) algorithm by incorporating the infinity norm of regression-based summary statistics to support the analysis of both binary and continuous outcomes for whole-genome and whole-exome sequencing studies. Compared to scan-based algorithms and the knockoffScreen method, dBiRS demonstrates greater power while maintaining strict control over family-wise error rates (FWER) and false discovery rates (FDR). Furthermore, dBiRS enables parallel computing of block-wise results, which are then aggregated through a central machine to ensure both detection accuracy and computational efficiency. Empirical studies further demonstrate its robustness and adaptability, even under conditions of signal decay.

Further extensions of dBiR could explore the integration of functional annotations. While the current framework effectively detects signal regions without annotations, their incorporation may provide additional biological context and improve detection accuracy for rare or weak signals in less characterized regions of the genome. Our framework is also flexible to incorporate various functional annotations by directly adding weights into SNP-level summary statistics, where weights can be estimated using annotation principal components (Zhou et al., 2023), CADD (Kircher et al., 2014) or other methods. Additionally, further validation of the identified associations through experimental studies is necessary to elucidate their functional relevance and causal mechanisms.

Detailed theoretical results on FWER control and detection accuracy are presented in the Appendix, while the proofs of the theorems are provided in the Supplementary Material.