Let a set of n SNPs have p-values P _i , i = 1, … , n. The TFisher statistic tests the genetic association between a trait and the SNP set by combining these p-values while allowing for a general scheme of truncation and weighting (Zhang et al., 2020b): T n τ 1 , τ 2 = − 2 ∑ i = 1 n log P i τ 2 I P i ≤ τ 1 , (1) where I () is the indicator function, τ ₁ > 0 is a truncation parameter that includes p-values equal or smaller than τ ₁ into the statistic, and τ ₂ > 0 is a weighting parameter for selected p-values. When τ ₁ = τ ₂ = 1, the TFisher statistic combines all p-values, which is the classic Fisher’s combination statistic: T n ( 1,1 ) = − 2 ∑ i = 1 n ⁡ log P i . The TFisher family include the statistic of truncation-product method (TPM) (Zaykin et al., 2002, 2007): T n ( τ 1 , 1 ) = − 2 ∑ i = 1 n ⁡ log P i I ( P i ≤ τ 1 ) , i.e., a special case of the TFisher with τ ₂ = 1. Our previous study has shown that statistical power and computation efficiency can be improved by weighting the truncated p-values through τ ₂. An optimality can be reached at τ ₁ = τ ₂ = τ ∈ (0, 1], which gives the soft-thresholding statistic: T n s τ = 2 ∑ i = 1 n max − log P i τ , 0 . An analogous version of the TPM is the rank-truncation product (RTP) method (Dudbridge and Koeleman, 2003). Let P ₍₁₎ ≤ ⋯ ≤ P _(n) be the ordered input p-values. The RTP statistic is R T P = − 2 ∑ i = 1 k ⁡ log P ( i ) for some predetermined k. The RTP statistic can also be written in consistency with the TPM with τ ₁ = P _(k). Calculating the p-value of the RTP is more challenging, especially for SNP p-values are dependent due to the LD.

The TFisher is a flexible framework to maximize the detection of SNP-set associations over a broad spectrum of signal patterns. Different signal patterns are in favor of different truncating and weighting parameters. For example, when association signals are dense, more SNP p-values should be included in the test statistic by large τ ₁ and τ ₂ so that the test is closer to Fisher’s method. Dense signals happen under the polygenic model with a substantial number of causal SNPs, or when the LD is strong so that many SNPs in LD with the causal SNPs also show association signals. On the other hand, if association signals are sparse (i.e., only a small number of SNP p-values are linked to the causal genetic factor), the smallest SNP p-values should be included in the statistic by small τ ₁ and τ ₂.

In reality, the signal patterns are often unknown and substantially differ over traits and loci. Therefore, we rely on the data-adaptive omnibus testing procedure to automatically select appropriate parameters. Specifically, we consider a discrete search domain over {(τ _1k , τ _2k ), k = 1, … , K}, where K is the total number of τ values to search on. Denote P(k) the test p-value of T _n (τ _1k , τ _2k ). The omnibus statistic is defined as the smallest P(k), which indicates the maximal association evidence for the whole SNP-set: oTFisher_minp = min k = 1 , … , K P k . (2) Moreover, we define a second omnibus test by Cauchy combination test (CCT) of P(k)’s (Liu and Xie, 2020): oTFisher_cct = 1 K ∑ k = 1 K tan 0.5 − P k π . (3) The summands of oTFisher_cct are the transformation of P(k)’s by the inverse cumulative distribution function (CDF) of the standard Cauchy distribution. Because of the heavy tail of Cauchy distribution, oTFisher_cct is dominated by p-values closer to 0 or 1. In practice, we truncate P(k) = 1 to be 0.9 so that oTFisher_cct is dominated by small p-values and performs similarly to oTFisher_minp. Note that since P(k) depends on the LD (see its calculation below), the oTFisher implicitly accounts for the LD information.

We can utilize the oTFisher as a procedure to screen for important SNPs. The oTFisher procedure has three steps:

1) SNP-set testing: Identify the significantly associated SNP-sets by their oTFisher p-values ≤ α / g , where α is the adjusted significance level, g is the number of SNP sets (e.g., genes) studied simultaneously.

A natural choice of the threshold is p ^⋆ = τ*, where τ * ≡ τ 1 k * corresponding to the oTFisher_minp in (2) (i.e., P (k*) is the minimal P(k)). Meanwhile, P(k) could have similar values over different k. To be conservative and reduce the false discoveries, we recommend p ^⋆ = min{τ*, 0.1} (denoted by oTFisher_r as a restricted version).

In practice, SNP screening is commonly based on the Bonferroni procedure or the Benjamini–Hochberg (BH) procedure:

• Bonferroni procedure: The screened SNPs are those with their p-values less than p ^⋆ = α/L, where L is the total number of SNPs.

The SNPs screened by Bonferroni and BH are validated in the same way as the validation stage for the oTFisher.

There are two potential benefits of utilizing the oTFisher procedure over Bonferroni and FDR procedures. First, as a set-testing method, the oTFisher can potentially increase the discovery of weakly associated SNPs. It is because the SNP-set test can detect the collective existence of weak genetic effects that are indistinguishable from individual SNPs (Donoho and Jin, 2004; Wu et al., 2014; Jin and Ke, 2016). Therefore, the oTFisher could better reveal SNPs with weak genetic effects than the Bonferroni and FDR procedures, which only rely on individual SNP tests. Second, τ* is influenced by the proportion of genetic signals (Zhang et al., 2020b). This information could also contribute to identifying important SNPs.

Following the literature (Brown, 1975; Zhang and Wu, 2022b), we account for the dependence of SNP p-values by assuming that the vector of their z-score statistics Z = ( Z 1 , … , Z n ) ′ is approximately normal: Z ∼ N μ , Σ , (4) where the mean vector μ corresponds to the association hypotheses: H ₀: μ = 0, i.e., no SNPs are associated, and H ₁: μ ≠ 0, i.e., at least one SNP is associated. The correlation matrix Σ is assumed to be estimable but otherwise arbitrary. These assumptions are reasonably satisfied in practice when the sample size is reasonably large (e.g., by the linear model-based association tests (Shao, 2010)). As one example, the estimation of Σ among the marginal score statistics is given in Section 3.1. For analyzing GWAS summary data where the individual-level genotype data are unavailable, Σ can often be estimated by the LD matrix based on reference genome panel data, such as the 1,000 Genome and the UK10K projects (Hu et al., 2013). Although most GWAS summary data contains two-sided p-values, we allow they are one-sided for the completeness of statistical development:  One - sided:  P i = Φ ̄ Z i ;  Two - sided:  P i = 2 Φ ̄ | Z i | , (5) where Φ ̄ ( x ) = P ( N ( 0,1 ) > x ) denotes the survival function of N (0, 1).

At given τ ₁, τ ₂ and n, the TFisher statistic T _n (τ ₁, τ ₂) in Eq. 1 has a point probability mass at 0: p 0 = P ( T n ( τ 1 , τ 2 ) = 0 ) = P ( min i P i > τ 1 ) corresponding to all p-values are truncated. Define s = −2 log (τ ₁/τ ₂) ≥ 0 when τ ₂ ≥ τ ₁, and s = −2n log (τ ₁/τ ₂) < 0 when τ ₂ < τ ₁. In either case T _n (τ ₁, τ ₂) ≥ s and its distribution is a mixture of point mass at 0 and a continuous distribution defined in [s, ∞). That is, T _n (τ ₁, τ ₂) ∼ p ₀ ⋅ 0 + (1 − p ₀) ⋅ T′, where T′ denotes an appropriate continuous random variable.

The exact value of p ₀ is easy to calculate under normality in Eq. 4. The exact distribution of T′ is challenging to obtain. We propose to use the gamma distribution model to approximate it for a few reasons. First, the model is consistent with the distribution of the TFisher under independence, which is a weighted gamma distribution (Zhang et al., 2020b). Second, the literature has been using gamma distribution to approximate Fisher’s method under dependence (Brown, 1975; Zhang and Wu, 2022b). Third, when the shape parameter of the gamma distribution is large, it converges to the normal distribution, which is appropriate for the TFisher statistic when n is large (see details below). Overall, gamma distribution provides a flexible and straightforward distribution model, vital for computational speed and accuracy.

Specifically, we approximate the distribution of the TFisher statistic by T n τ 1 , τ 2 ≈ D X ∼ p 0 ⋅ 0 + 1 − p 0 ⋅ X ′ , (6) where X′ − b ∼Γ(a, θ), the gamma distribution with shape parameter a and scale parameter θ. We consider a shift parameter b so that X′ ∈ [b, ∞). Let F _{Γ(a,θ) (} x) denote the CDF of Γ(a, θ). The CDF of X is F X x = P X ≤ x = p 0 I x ≥ 0 + 1 − p 0 F Γ a , θ x − b . (7) Based on (7), the p-value of the TFisher at an observed statistic t is P T n τ 1 , τ 2 > t ≈ 1 − F X t . (8) We discuss the methods to calculate F _X (t) in the following.

Method 1: The generalized Brown’s method (GB): This method follows the essential idea of Brown’s method (Brown, 1975) to match the first two moments of T and X. Specifically, we set the shift parameter b = s so that X′ and T′ have the same domain. The parameters a and θ are determined by matching the means and the variances of T _n (τ ₁, τ ₂) and X. Denote μ _T = E (T _n (τ ₁, τ ₂)) and σ T 2 = V a r ( T n ( τ 1 , τ 2 ) ) . We have μ T = μ X = 1 − p 0 a θ + b σ T 2 = σ X 2 = 1 − p 0 a θ 2 + p 0 a θ + b 2 ⇒ a = μ T − b 1 − p 0 2 1 − p 0 σ T 2 − p 0 μ T 2 θ = 1 − p 0 σ T 2 − p 0 μ T 2 1 − p 0 μ T − b 1 − p 0 . Note that the gamma approximation is consistent with the asymptotic normal distribution of T _n (τ ₁, τ ₂) for large n by the Central Limit Theorem (CLT). Specifically, when n → ∞, p ₀ → 0, so μ _T = μ _X ≈ aθ + b and σ T 2 = σ X 2 ≈ a θ 2 . Because Γ(a, θ) ≈ N (aθ, aθ ²) for large a, the distribution model in (6) leads to T n ( τ 1 , τ 2 ) ≈ D Γ ( a , θ ) + b ≈ N ( a θ + b , a θ 2 ) ≈ N ( μ T , σ T 2 ) . However, for finite n, the distribution model in Eq. 6 is more accurate for p-value calculation.

Straightforward calculation gives μ _T = 2nτ ₁ (1 − log  τ ₁ + log  τ ₂). For the variance σ T 2 , we deduce its analytical formula given in Lemma 1 in Supplementary Material. The formula involves a summation of infinite terms. However, in practice, a summation of two or three terms over k would give sufficient accuracy for σ T 2 (Zhang and Wu, 2022b). The proof is based on Mehler’s theorem (Patel and Read, 1996) and is given in Supplementary Material.

Method 2: Skewness-kurtosis-ratio matching method (SKRM). Accurate calculation of small p-value highly depends on the precise approximation of the right tail of the null distribution. In this method, we do not require the shifting parameter b = s but treat it as additional freedom to capture the right-tail information of the TFisher statistic. That is, in addition to the first two moments, we further match the skewness-kurtosis ratios of T _n (τ ₁, τ ₂) and X. Engaging higher-order moments could provide more flexibility in the distribution and thus improve the accuracy of p-value calculation. In particular, matching the skewness-kurtosis ratios is a cost-efficient method–it captures two higher moments using only one extra parameter.

Specifically, let γ T ′ = E ( T ′ − μ T ′ ) 3 / σ T ′ 3 and κ T ′ = E ( T ′ − μ T ′ ) 4 / σ T ′ 4 be the skewness and kurtosis of T′, and γ _X′ and κ _X′ be the skewness and kurtosis of X′, respectively. By matching the ratio between skewness and excess kurtosis γ T ′ κ T ′ − 3 = γ X ′ κ X ′ − 3 , we can obtain a simple closed form a = 9 γ T ′ 2 ( κ T ′ − 3 ) 2 . Subsequently, by matching the mean and variance μ _T′ = μ _X′ and σ T ′ 2 = σ X ′ 2 , we have θ = σ T ′ a and b = μ T ′ − σ T ′ 2 a . After a, θ and b are determined, the p-value of the TFisher can be calculated by Eqs 7, 8.

Exact calculation of γ _T′ and κ _T′ would be intricate due to the complexity of the high moments of the summational terms in the TFisher statistic. We rely on simulation by Eq. 4 to obtain these values. The number of simulations needed for estimating these parameters is much smaller than that required for obtaining a small empirical p-value directly. Therefore, the SKRM is still computationally more efficient than getting p-values solely based on the re-sampling strategy.

Method 3: Hybrid method. To balance computational speed and accuracy, we recommend a simple two-stage calculation of TFisher’s p-value. Because the GB is fast and accurate for controlling the type I error rate at α ≥ 0.01 (see the numerical results below), the GB is applied in the first stage. If the obtained p-value is less than 0.01 (the chance is about 1% under the null), the SKRM method will obtain the final p-value in the second stage. With a single core of 2.80 GHz AMD EPYC 7543 CPU and 20G memory, the computation times for calculating TFisher’s p-values for SNP sets of 30/50/100/200 SNPs are about 0.07/0.13/0.30/0.74 s by the GB method (implemented in R version 4.2.0). Correspondingly, the expected times by the hybrid method (assuming 1% chance of engaging the SKRM method that takes 10⁵ simulations to obtain γ _T′ and κ _T′ values) are 0.09/0.15/0.33/0.78 s.

For the oTFisher_minp in Eq. 2, we can apply asymptotic distribution to approximate its p-value by 1 − Φ R Φ ̄ − 1 p 0 , … , Φ ̄ − 1 p 0 , (9) where p ₀ is the observed statistic of oTFisher_minp, Φ _R denotes the CDF of a multivariate normal distribution with mean zero and correlation matrix R. We obtain R by scaling Ω, the covariance matrix of T _n (τ _1k , τ _2k ), k = 1, … , K, given by Lemma 2 in Supplementary Material. That is, R = ΛΩΛ with the diagonal matrix Λ = diag { 1 / Ω k k , k = 1 , … , N } .

As for oTFisher_cct in Eq. 3, following the property of the CCT, its distribution is robust to the correlations as long as T _n (τ _1k , τ _2k )’s are roughly normal distributed. This requirement is justifiable because T _n (τ _1k , τ _2k ) is in the format of summation so that it is roughly normal distributed by the CLT when n is moderately large and τ ₁ is not too small. Denote the observed statistic by cct _o . We directly apply the result by (Liu and Xie, 2020) to approximate its test p-value by P ( o T F i s h e r _ c c t > c c t o | H 0 ) ≈ 1 2 − tan − 1 ( c c t o ) / π . For the computation time of obtaining the oTFisher’s p-value, the dominant part is to get the p-values of the TFisher statistics T _n (τ _1k , τ _2k ), k = 1, … , K.

Simulations were applied to verify the accuracy of the p-value calculation, statistical power, and SNP-screening performance of the oTFisher procedures. The genotype data were generated by the Cosi2 package (Shlyakhter et al., 2014). Specifically, 1,290 haplotypes were generated according to a coalescent model based on chromosome 1 of the European population. Two haplotypes were randomly picked with replacement to form the genotypes of one diploid individual. In each simulation, we obtained SNPs of N individuals. Both rare variants (0.05% < MAF < 5%) and common variants (MAF ≥5%) were considered.

We simulated continuous and binary traits by the regression and the logit model, respectively: Y k = G k ⋅ ′ β + X k ⋅ ′ γ + ϵ k ,  where  ϵ k ∼ i . i . d . N 0,1 , logit P Y k = 1 = G k ⋅ ′ β + X k ⋅ ′ γ , (10) where Y _k quantifies the phenotypic trait of the kth subject, k = 1, … , N, with the sample size N. G k ⋅ = ( G k 1 , … , G k n ) ′ is the genotype vector of n SNPs, X k ⋅ = ( X k 1 , … , X k l ) ′ is the vector of l controlling covariates. The nonzero elements of β = ( β 1 , … , β n ) ′ are the causal genetic effects of the corresponding SNPs. The SNP-set analysis concerns testing the global hypotheses H 0 : β = 0  versus  H 1 : β ≠ 0 . We mimic a balanced case-control study for the binary traits. That is, a large number of outcomes were generated based on the probability of the logit model, then we randomly selected subjects so that the numbers of cases and controls are N/2 each.

Based on the simulated data, we calculate the marginal score test statistic following literature (McCullagh and Nelder, 1989; Schaid et al., 2002; Barnett et al., 2017). Specifically, the score of the ith SNP is M i = ∑ k = 1 N G k i ( Y k − Y ̃ k ) , i = 1 , … , n , where Y k ̃ are the fitted trait values by the maximum likelihood estimation under H ₀. It can be shown that under H ₀, M = ( M 1 , … , M n ) ′ → D N ( 0 , Σ ) , as N → ∞. The covariance matrix Σ can be estimated by Σ ̂ = G ′ W G − G ′ W X ( X ′ W X ) − 1 X ′ W G , where (G, X) is the design matrix corresponding to (10), and W is a diagonal matrix: W = σ ̂ 2 I for continuous trait (where σ ̂ 2 is the estimate of the residual variance); W = diag { Y ̃ k ( 1 − Y ̃ k ) , k = 1 , … , N } for binary trait. Each M _i is standardized to get the marginal score statistic Z i = M i / Σ ̂ i i → D N ( 0,1 ) as N → ∞ under H ₀. The correlation matrix of Z = ( Z 1 , … , Z n ) ′ is estimated by D ̂ Σ ̂ D ̂ , where D ̂ = diag { Σ ̂ i i , i = 1 , … , n } . Thus, the marginal score statistics satisfy the assumption (4) asymptotically. We used the two-sided SNP p-values in Eq. 5. For the rare-variant analysis of binary traits, the saddle point approximation (SPA) was applied to obtain SNP p-values, which corrects the bias due to the unbalanced distribution of rare variants’ genotype data (Dey et al., 2017).

Under the null, quantitative and binary trait values were generated by setting β = 0 and Z = 1 in Eq. 10. We simulated 10⁷ oTFisher statistics with the adapting domain τ 1 = τ 2 = τ ∈ T = { 0.001 , 0.005 , 0.01 , 0.05 , 0.1 , 0.2 , 0.5 , 0.7 , 1 } . The simulations included 100 randomly generated genotype data to mimic various minor allele frequencies and LD structures of SNPs in various genes. The empirical type I error rate was obtained by proportionating all calculated p-values smaller than a given nominal level α. Table 1 lists the ratios of empirical type I error rates and the nominal α levels under common and rare variants for quantitative and binary traits. A ratio around 1 indicates accurate calculation. A more stringent α is harder to control. The GB method well controls type I error rate up to α ≥ 0.005, but becomes liberal at smaller α, where the SKRM method (using 10⁵ simulations to obtain the third and fourth moments) controls the error much better. Therefore, the hybrid method combining the GB and the SKRM balances accuracy and computational speed. Consistent simulation results with various n and T settings (with starting τ = 10^–5, 10^–3 or 10^–2 in T ) are given in Supplementary Tables S1–S4 in Supplementary Material.

Through simulations, we assessed the statistical power of the oTFisher tests in comparison with other p-value combination tests in GWAS summary data analysis: the GATES (extended Simes procedure, using GATES2 function in R library aSPU (Li et al., 2011)), the CCT (i.e., the ACAT with equal weights), Fisher’s method, the soft-thresholding TFisher at fixed τ ₁ = τ ₂ = 0.05, the adaptive TPM (ATPM), and the adaptive RTP (ARTP). The ATPM follows (1) and (2) with fixed τ ₂ = 1, and adapts over τ 1 ∈ T . The ARTP adapts over τ ₁ ∈ {P _(k)}, where k ∈ n T (rounding to the nearest integers) to be consistent with the oTFisher and the ATPM.

We considered that the causal SNPs were randomly located, and their effect sizes were given by the nonzero elements of β in Eq. 10. For a fair comparison, we empirically controlled the type I error rate α to avoid potentially unavailable or inaccurate p-value calculation for some tests. For example, there are no p-value calculation methods for the RTP and the ARTP under dependence. We got the critical value of one statistic by the upper 100α% percentile of its values generated from 10,000 simulations under the null. The statistical power was obtained by the percentage surpassing the critical value among the statistics generated from 1,000 simulations of the alternative. For the RTP and ARTP, we applied a one-level simulation algorithm consistent with literature (Yu et al., 2009) except that we directly simulated the Z-scores for faster computation instead of permuting the genotype data (details see Supplementary Material).

Statistical power was systematically studied under various settings regarding trait type (quantitative or binary), SNP type (common or rare), the number of causal SNPs m, genetic effect size β, SNP-set size n, sample size N, and type I error rate α. Figure 1 shows the power comparison under binary traits from common SNPs at α = 0.005 with n = 100 and N = 10,000. A few interesting observations can be made. First, the GATES and the CCT have similar performances. They are advantageous when causal SNPs are sparse (e.g., m ≤ 3) and their effects are strong. Fisher’s method shows an opposite pattern—it is preferred if causal SNPs are dense, especially when effects are weak. These patterns are consistent with literature results (Zhang and Wu, 2022a). TFisher at fixed τ ₁ = τ ₂ = 0.05 is more robust over sparse and dense causal SNPs. However, it could still be less satisfactory (e.g., when m ≤ 3). The oTFisher_cct and oTFisher_minp are similar; their power is the best in most scenarios, showing an overall advantage over unknown genetic architectures. In comparison with the ATPM and the ARTP, Figure 2 shows that the oTFisher is uniformly better than the ATPM. This observation is well-supported by a theoretical optimality study (Zhang et al., 2020b). The ARTP and the oTFisher have very similar power. Meanwhile, our p-value calculation algorithm provides a practical advantage for applying the oTFisher over the ARTP in computation. These comparison patterns remain similar for quantitative traits, rare variants, and different α levels. Comparisons under other settings are given in Supplementary Figures S1–S13 in Supplementary Material.

We studied the performance of SNP-screening procedures measured by the accuracy of detecting causal SNPs. Because non-causal SNPs in LD with causal SNPs also show statistical associations, the study focused on rare-variant analysis with weak LDs for simplicity. To mimic a gene-based SNP-set analysis, we simulated L = 1,000 SNPs (with the LD r ² < 0.3) in g = 10 genes of equal size. Two causal genes contained causal SNPs with random locations. The continuous and binary traits were obtained using models in Eq. 10 that included all causal SNPs. We systematically varied the genetic effect β and the proportion of causal SNPs in the two causal genes. The sample size N = 1,000; the cases and controls were balanced for binary traits.

We considered accuracy by the sensitivity, specificity, and balanced accuracy (BA, the average of sensitivity and specificity) based on the true positives (TP, the picked SNPs that are causal), false positives (FP, the picked SNPs that are non-causal), true negatives (TN, the unpicked SNPs that are non-causal), and false negatives (FN, the unpicked SNPs that are causal). These numbers are determined after defining the “picked” and “unpicked” SNPs. At the screening stage of the oTFisher procedure described in Section 2.2, we consider the s candidate SNPs as being picked from in total t SNPs in the significant genes; the rest t − s SNPs are unpicked. At the validation stage, the s ₁ validated SNPs are picked, and the rest t − s ₁ SNPs in the significant genes are unpicked. The accuracy measures were averaged over 1,000 simulations.

The oTFisher procedure was compared with Bonferroni and BH procedures for SNP screening. We further considered an oracle procedure:

• Oracle procedure: Assume the number of causal SNPs m is known, the oracle, i.e., the best possible, SNP screening procedure is to pick SNPs by setting p ^⋆ = P _(m), the mth smallest SNP p-value.

Certainly, m is unknown in reality, so this procedure is a hypothetically optimal procedure serving as an indicative accuracy for comparison purposes. The validation process is the same for all procedures.

Figure 3 compares the BA of the screened SNPs and validation of them. Overall, in both screening and validation stages, the oTFisher_r gave a higher BA than Bonferroni and the BH, sometimes even the Oracle. In the screening stage, the oTFisher using p* = τ* had higher sensitivity but lower specificity. Restricting the threshold to p ^⋆ = min{τ*, 0.1} significantly increased the specificity. The validation stage helped further control the type I error. For comparisons of sensitivity and specificity and results under more settings, see Supplementary Figures S14–S21 in Supplementary Material.

We studied the genetic associations between genes and traits using the SNP sets grouped by genes. Four methods were applied. First, the Bonferroni procedure represents the single-SNP method applied in the original GEFOS studies. The top-hit genes contained significant SNPs with p-values less than Bonferroni adjusted threshold: 0.05 divided by the total number of SNPs in all genes. Second, we applied two broadly applied SNP-set methods using SNP p-values: the ACAT (Liu et al., 2019) and the GATES. Third, we applied our proposed oTFisher method with the adapting domain τ 1 = τ 2 = τ ∈ T = { 1 0 − 5 , 0.001 , 0.005 , 0.01 , 0.05 , 0.1 , 0.2 , 0.5 , 0.7 , 1 } .

The oTFisher_cct and the oTFisher_minp yielded similar results, so the former is reported below for simplicity. The Q-Q plots of the gene-based association p-values are given in Supplementary Figure S24. The top-hit genes are given in the supplementary file top_gene.xlsx. The overlaps of the top-hit genes and the top-hit literature genes among these four analysis methods are given in Supplementary Figures S25–S26.

Figure 4 summarizes the numbers of top gene hits, from which we can make a few interesting observations. First, the Bonferroni procedure systematically led to fewer top hits and literature genes than the other methods. The result indicates that SNP-set tests could have higher statistical power than single-SNP analysis in detecting disease genes. Second, the oTFisher yielded similar or more gene hits, and most of the hits are literature genes indicating a reliable discovery and potentially higher statistical power. In particular, the oTFisher could have the advantage of detecting polygenic genes that contain relatively dense genetic signals. For examples, the oTFisher detected significantly more genes in the studies of GEFOS2017_TBBMD and GEFOS2020_FALLS, where top-hit genes often contain multiple SNPs with relatively small p-values (see Supplementary Tables S9 and S10 for the distribution of SNP p-values within these top-hit genes). The polygenic genetic architecture is possible for complex human diseases, including the BMD-related traits (Kemp et al., 2017; Morris et al., 2019).

We carried out gene-set enrichment analyses (GSEA) for the top-hit genes identified by oTFisher. Based on gene ontology (GO) and KEGG pathways, the analysis was conducted by KEGG Orthology-Based Annotation System intelligent version (KOBAS-i) (Bu et al., 2021). The analysis identified significant GO terms and biological pathways enriched in the top-hit genes at the corrected significance level of 0.05. The GEFOS data-specific results are summarized in supplementary files top_GOs_pathways_study-specific.xlsx. These GO terms and biological pathways are often related to the bone and skeletal system and are consistent with the osteoporosis pathways reported in literature (Guo et al., 2019). Furthermore, we considered the BMD as a general trait and carried out the GSEA by pooling 173 top-hit genes from eight studies (except for falling risk). Twenty-one significantly enriched pathways were obtained and clustered into three networks based on their correlations by the cirFunMap plot (Bu et al., 2021). Figure 5 visualizes the clusters (ranked by the median p-value of the enriched pathways within each cluster). The first cluster gives a major network containing 12 pathways: Wnt signaling pathway (hsa04310), breast cancer (hsa05224), hepatocellular carcinoma (hsa05225), pathways in cancer (hsa05200), gastric cancer (hsa05226), basal cell carcinoma (hsa05217), signaling pathways regulating pluripotency of stem cells (hsa04550), proteoglycans in cancer (hsa05205), hippo signaling pathway (hsa04390), human papillomavirus infection (hsa05165), Cushing syndrome (hsa04934), and mTOR signaling pathway (hsa04150). The top two significant pathways, the Wnt signaling and breast cancer pathways, have been reported in literature (Guo et al., 2019). The large cluster here provides a networking context for them. The second cluster, containing pancreatic cancer (hsa05212) and colorectal cancer (hsa05210), is also connected with the first cluster. The third cluster, containing the prolactin signaling pathway (hsa04917) and rheumatoid arthritis (hsa05323), is independent of the rest. Details of the significant pathways and their clusters for the BMD traits are given in supplementary files top_pathways_BMDs.xlsx.

The top-hit novel genes included in the enriched GO terms and pathways are likely disease genes that influence the corresponding functionalities. In particular, we obtained three top-hit novel genes in the 21 significantly enriched pathways obtained by pooling 173 top-hit genes of the BMD traits. Gene HSPG2 (chr1: 22148724–22263790, oTFisher p-value 1.18E-07) is included in a significantly enriched pathway of proteoglycans in cancer (hsa05205, corrected enrichment p-value 0.0047). It was shown to be associated with segregating developmental dysplasia of the hip (Basit et al., 2017). Gene MAP3K12 (chr12: 53874275–53893444, oTFisher p-value 6.49E-09) is included in the significantly enriched MAPK signaling pathway (hsa04010, corrected enrichment p-value 0.0048). It is related to lissencephaly type 3 - metacarpal bone dysplasia and infantile osteopetrosis with neuroaxonal dysplasia in the Open Targets Genetics (Ghoussaini et al., 2021). Gene PRKAG1 (chr12: 49396054–49412629, oTFisher p-value 2.22E-06) is included in the significantly enriched Apelin signaling pathway (hsa04371, corrected enrichment p-value 0.0049). It is related to bone marrow failure syndrome in the Open Targets Genetics. Novel top-hit genes contained in the enriched GOs and pathways from GEFOS data-specific results are summarized in the supplementary file novel_genes_in_top_GOs_pathways_study-specific.xlsx. More discussion of them is given in Supplementary Material, including genes connected to relevant traits such as osteoarthritis, osteosarcoma, and bone metastasis.

Gene-based analysis has the limitation of a small coverage of the genome. For a whole-genome association study, we grouped and analyzed SNPs by haplotype blocks (haploblock estimation by PLINK (Chang et al., 2015) is detailed in Supplementary Material). The Q-Q plots and the Manhattan plots are given in Supplementary Figures S27 and S28 in Supplementary Material. Overall, genomic inflation is reasonably controlled. For various SNP-set methods, Figure 6 shows the number of top-hit blocks. Compared with other methods, the oTFisher generated more top-hit blocks and novel blocks (i.e., top-hit blocks that do not overlap literature genes or SNPs). More details on the top-hit blocks and their corresponding SNPs and genes are summarized in the supplementary file top_haploblocks.xlsx.

The haploblock-based analysis provided complementary results to the single-SNP analysis and gene-based analysis. As a SNP-set analysis method, haploblock analysis could detect additional disease SNPs over single-SNP analysis. For example, in the study GEFOS2012_FN, the top-hit blocks by the oTFisher discovered 56 literature SNPs that the single-SNP analysis failed to detect (since their p-values do not pass the genome-wide significance level). The haploblock analysis could also map additional disease genes over the gene-based analysis. Supplementary Table S11 summarizes the numbers of top-hit blocks that can map to literature or novel genes. According to the results, the haploblock analysis found some literature genes that were not among the top hits of the gene-based analysis. For example, in the study GEFOS2012_FN, six literature genes were mapped by top-hit haploblocks but were not discovered by the gene-based analysis: ATXN7L3, AXIN1, CPED1, FUBP3, LOC100272217, and SOX6. For all GEFOS studies, we had 27 literature genes and 119 literature SNPs (63 of them are outside gene regions) detected by haploblock analysis but not by gene-based analysis. Supplementary Figure S29 lists the numbers of literature genes and SNPs found by haploblock analysis versus gene-based analysis. Furthermore, the top-hit blocks (including single-SNP blocks) contained all GEFOS-reported significant SNPs, indicating no information lost compared to the original GEFOS studies.

The top-hit blocks from eight BMD studies contain 286 novel blocks; 255 of them are outside of genes (detailed information is given in the supplementary file novel_haploblocks.xlsx). By epigenetic annotation (Haploreg v3 (Ward and Kellis, 2016)), 58 of the novel blocks (representing 43 non-overlapping loci) co-locate with strong enhancers of literature genes. Therefore, these novel blocks are of interest due to their functional connections. For example, a top-hit block chr2:54643778-54645650 (oTFisher p-value 1.67E-16 by GEFOS2012_LS data) contains three SNPs rs13393949 (p-value 2.12E-06), rs4671215 (p-value 9.36E-08), and rs7560205 (p-value 1.72E-05) that locate at a strong enhancers of gene SPTBN1 in cells Huvec (umbilical vein endothelial cells) and NHEK (epidermal keratinocytes). Gene SPTBN1 is shown associated with heel bone mineral density (Kemp et al., 2017). The enhancer-located novel blocks are given in the supplementary file novel_haploblocks_enhancers_literaturegene.xlsx, and an extensive discussion of the related SNPs is given in Supplementary Material.

Real-data analysis results show that the SNPs-screening procedure by the oTFisher_r could likely yield more disease SNPs than the Bonferroni procedure. Specifically, we used GEFOS2012 data sets to collect screened SNPs and validate them by a large data of osteoporosis from the UK Biobank (15,133 cases and 426,942 controls at ages of 38–73 years) (Sudlow et al., 2015). Figure 7 shows the Venn diagrams of the screened SNPs by the oTFisher_r and the Bonferroni procedures in the haploblock-based analyses, which are also compared with the SNPs reported in the original GEFOS studies and the literature SNPs. For consistent comparison, the validated SNPs were defined by the significance level of 0.05 over the total number of unique SNPs from the screening step and the literature. As expected, all GEFOS-reported SNPs were contained by the sets of literature SNPs as well as the screened SNPs. The oTFisher_r replicated more literature SNPs than the Bonferroni both before and after the validation stage. For example, with GEFOS2012_FN data, all screened SNPs by the Bonferroni were included in the set of SNPs by the oTFisher_r, while the oTFisher_r screened 30 additional literature SNPs (among which five were validated). Therefore, the oTFisher_r likely has a higher chance of finding disease SNPs. Furthermore, over 30% of the screened SNPs by oTFisher_r were verified. The high validation percentage (compared to the expected percentage of no more than 5% under the null) indicates that the set of screened SNPs by oTFisher_r likely contains enriched disease SNPs. Consistent results for gene-based analysis are given in Supplementary Figures S30 and S31 in Supplementary Material.