Weighted multiple testing procedures in genome-wide association studies

Benjamini and Hochberg

The linear step-up procedure (BH) proposed by Benjamini & Hochberg (1995) to control the FDR at level α consists in rejecting all k null hypotheses corresponding to the k smallest p-values where $k=max\left(i\ge 0:p_{\left(i\right)}\le \frac{i\alpha }{m}\right)$, p_(i) being the ordered p-values. It has been shown that when the test statistics are PRDS (positive regression dependent on subset of null statistics), the BH procedure controls the FDR at level π₀α (i.e., FDR ≤ π₀α) (Benjamini & Hochberg, 1995; Benjamini & Yekutieli, 2001).

Qvalue

To improve the power of the BH procedure, several adaptive procedures have been introduced in which the proportion of true null hypotheses π₀ is estimated from the data (Storey, 2002; Dalmasso, Broët & Moreau, 2005; Benjamini, Krieger & Yekutieli, 2006; Liang & Nettleton, 2012). One of the most used is the Qvalue procedure, in which a cubic spline based method is used to estimate the quantity $li{m}_{\lambda \to 1}{\hat{\pi }}_0\left(\lambda \right)$ where ${\hat{\pi }}_0\left(\lambda \right)=\frac{\#\left\{p_i>\lambda \right\}}{m\left(1-\lambda \right)}$ (Storey, 2002; Storey & Tibshirani, 2003).

Weighted Benjamini and Hochberg

The weighted BH procedure (wBH) was introduced by Genovese, Roeder & Wasserman (2006). It consists in assigning to each null hypothesis H_0,i a non-negative weight such that $\stackrel{m}{\sum _{i=1}}{w}_i=m$. Then, the BH procedure is applied by replacing p_i by $\frac{p_i}{w_i}$. This procedure has been proven to control the FDR.

False discovery rate regression

The false discovery rate regression (FDRreg) procedure introduced by Scott et al. (2015) is an adaptive procedure in which the proportion of true null hypotheses π₀ is estimated. However, this quantity is made dependent on the covariate, leading to: π₀(x_i) = Pr(H_i = 0|X_i = x_i) and $FDR\left(x_i\right)=\mathbb{E}\left(\frac{FP}{max\left(R,1\right)}|X_i=x_i\right)$ with π₀(x_i) representing weights specific to each hypothesis. Thus, noting z_i the test statistics, the two-components mixture model can be written: (2)${\displaystyle f\left(z_i\right)={\pi }_0\left(x_i\right)f_0\left(z_i\right)+\left(1-{\pi }_0\left(x_i\right)\right)f_1\left(z_i\right)}$

In this approach, the alternative density f₁(z_i) is taken to be a location mixture of null density which is assumed to be Gaussian. The mixing proportion of f₁(z_i) is fitted via a predictive recursion algorithm (Newton, 2002) and the model parameters in Eq. (2) are then estimated by an EM algorithm treating the mixing proportions of the alternative distribution as fixed. A fully Bayesian approach based on MCMC method is also proposed.

Science-wise false discovery rate

The science-wise false discovery rate (swfdr) procedure introduced by Boca & Leek (2018) is similar to FDRreg in the sense that π₀ and FDR are extended by conditioning on the covariate. However, while FDRreg jointly estimates π₀ and the FDR by assuming that the test statistics are normally distributed, swfdr first estimates the proportion of true null hypotheses, and then the FDR is obtained from a plug-in estimator. To estimate the proportions π₀(x_i), an approach similar to Qvalue is proposed but the ratio $\frac{\#\left\{p_i>\lambda \right\}}{m\left(1-\lambda \right)}$ is replaced by $\frac{\hat{\mathbb{E}}\left(1_{P_i>\lambda }|X_i=x_i\right)}{\left(1-\lambda \right)}$ where 𝔼(1_Pi>λ|X_i = x_i) is estimated from a logistic regression model.

Covariate adaptive multiple testing

The covariate adaptive multiple testing (CAMT) procedure introduced by Zhang & Chen (2020) is also based on the mixture model (Eq. (2)) with mixing proportions dependent on the covariate. However, this procedure relies on the local version of the FDR, the local false discovery rate (lfdr), introduced by Efron et al. (2001). The lfdr is defined as the posterior probability that a hypothesis is null given a specific pvalue: $lfdr\left(p_i\right)=Pr\left(H=0|P=p_i\right)=\frac{{\pi }_0{f}_0\left(p_i\right)}{f\left(p_i\right)}$.

From this definition, the FDR can be derived by the relationship FDR = 𝔼(lfdr|P ∈ Γ) where Γ is a rejection region for the pvalues (Efron et al., 2001; Dalmasso, Bar-Hen & Broët, 2007). In addition, the optimal decision rule can be written $lfdr\left(p_i\right)\le t\iff \frac{f_{1,i}\left(p_i\right)}{f_0\left(p_i\right)}\ge \frac{\left(1-t\right){\pi }_0\left(x_i\right)}{t\left(1-{\pi }_0\left(x_i\right)\right)}$.

The principle of CAMT is to replace the ratio $\frac{f_{1i}}{f_0}$ in the optimal decision rule by a surrogate function $h_i\left(p\right)=\left(1-k_i\right)p_i^{-k_i}$. Then, the parameters π₀(x_i) and k_i are estimated from an EM algorithm in order to find the optimal threshold t that allows the FDR to be controlled at the desired level.

Independent Hypothesis Weighting

The Independent Hypothesis Weighting procedure (IHW) was introduced by Ignatiadis et al. (2016). Here, the objective is to find optimal weights that maximize overall power. The basic idea is to divide hypotheses into G groups according to the ordered values of the covariate. Then, positive weights are assigned to each group g in order to maximize the number of rejections.

To avoid overfitting, the authors introduced a hypothesis splitting approach which consists in randomly splitting the m hypotheses into k folds independently of the pvalues and covariates. For each fold, an optimization problem is applied to the hypotheses of the k − 1 remaining folds in order to derive weights ${\tilde{\omega }}_g$ (g = 1, …, G) that maximize the overall power. Then, hypotheses of the held out fold lying in the group g are assigned weight ${\tilde{\omega }}_g$.

To make the optimization problem convex, the authors proposed using the Grenander estimator instead of the empirical cumulative distribution function. In addition, to solve the optimization problem, they added a regularization parameter λ such that ${\sum }_{g=2}^G\parallel w_g-w_{g-1}\parallel \le \lambda$ where λ > 0. This regularization parameter allows for weights of successive groups to be relatively similar.

Finally, once the weights are estimated, a standard wBH procedure is applied with the resulting weights vector.

Genotypes

Note G the genotype matrix with n lines corresponding to individuals (set to 2000) and m columns corresponding to SNPs (m ∈ {8000, 14000, 20000}). To code the genotype matrix, we considered an additive genetic model with 0 for a homozygous genotype for the reference allele, 1 for a heterozygous genotype, and 2 for a homozygous genotype for the alternative allele. Thus, G is a matrix of size n × m where G_ij ∈ {0, 1, 2}.

To mimic the linkage disequilibrium (LD) structure, we considered a model adapted from the work of Wu et al. (2009). The full genotype of each individual was generated from an m-dimensional multivariate normal distribution: $G_i^{\ast }\sim {\mathcal{N}}_m\left(0,\Sigma \right)$. The correlation matrix Σ is block-diagonal with blocks of size B = 10 and within each block, all variables are equicorrelated at level ρ. The different values considered for ρ are: 0, 0.10, 0.20, 0.35, 0.5 and 0.75.

To obtain the genotypes, these continuous variables were discretized following the Hardy Weinberg equation: p² + q² + 2pq = 1 where, for each SNP, p is the frequency of one of the two possible alleles and q = 1 − p.

From here, we arbitrarily set p as the MAF, so that p ≤ q. The MAF of the m₀ non-causal variants were generated from a uniform distribution between 0.01 and 0.5 (U[0.01, 0.5]). The m₁ causal variants were divided into four distinct subsets in which the MAF were generated from the following distributions:

• Group 1 (Rare SNPs): U[0.01, 0.05]

• Group 2 (Medium-Rare SNPs): U[0.05, 0.15]

• Group 3 (Medium SNPs): U[0.15, 0.25]

• Group 4 (Common SNPs): U[0.30, 0.40]

The number of SNPs in each subgroup of causal variants was obtained by the Euclidean division quotient of m₁ by 4 with m₁ ∈ {5, 10, 15, 20, 25, 50, 100, 150}. The remainder of this division was added to group 4. The effect size simulation is described in the next section.

In the end, we set for each SNP:

• G_ij = 2 if $G_{ij}^{\ast }<q_{p^2,N\left(0,1\right)}$,

• G_ij = 1 if $q_{p^2}<G_{ij}^{\ast }<q_{{\left(1-p\right)}^2,N\left(0,1\right)}$,

• G_ij = 0 if $q_{{\left(1-p\right)}^2,N\left(0,1\right)}<G_{ij}^{\ast }$,

where q_.,N(0,1) is the quantile function of the standard normal distribution.

Phenotypes

Let Y = (y₁, …, y_n) be the studied phenotype. In our simulation setting, we considered quantitative and binary variables corresponding to studies for quantitative trait and case-control design, respectively.

Larger numbers of tested hypotheses

To evaluate the different procedures with more realistic numbers of hypotheses in GWAS, we also simulated datasets with larger m and m₁ values (m ∈ {100000, 200000, 500000} and m₁ ∈ {100, 150, 250}) with R² = 0.5. However, due to excessive computational time for generating the data, we only considered the independent case with quantitative traits.

Weights and covariates

In this work, we sought to prioritize the detection of rare variants having strong genetic effects in GWAS. Thus, for wBH, we set the weights analogously to wBHa by considering $w_i=\frac{m}{{\sum }_{j=1}^m\frac{1}{x_j}}\times \frac{1}{x_i}$ (with x_i being the MAF). For the other weighting procedures, the MAF was used as the informative covariate.

However, as previously mentioned, other informative covariates may be considered. For example, we illustrate the use of the proposed method to prioritize common variants by replacing MAF by 1/MAF for all methods (except for CAMT for which we used log(1/MAF) to avoid computational problems in the EM-algorithm due to large values of the covariate). In addition, to evaluate whether the proposed method is robust when the covariate is completely uninformative, we applied all procedures with a covariate drawn from a uniform distribution between 0 and 1 (U[0, 1]).

Package versions

Packages, their versions, and functions used for all analyses are displayed in Table 3.

Overall power

To evaluate the ability of the procedures to detect true associations, for each configuration, we estimated the average power $\mathbb{E}\left(\frac{TP}{m_1}\right)$ by the empirical mean (and its corresponding standard error) of the numbers of true discoveries over the 500 simulated datasets divided by m₁.

However, for correlated statistics, the definitions of true and false positives are ambiguous since one single genetic marker influencing the phenotype may lead to multiple significant results for correlated markers (Benjamini, Krieger & Yekutieli, 2006; Siegmund, Yakir & Zhang, 2011; Brzyski et al., 2017). To handle this problem, we estimated the power (and the FDR) in correlated datasets by considering clusters of correlated SNPs as units of interest. The clusters were defined according to an estimated correlation coefficient threshold of 0.8.

Power in subgroups

To evaluate the performance of the methods specifically in the subgroup of rare variants, we also estimated the average power in each subgroup $\mathbb{E}\left(\frac{T{P}_g}{{m_1}_g}\right),g=\left(1,2,3,4\right)$ (where g = 1 corresponds to the subgroup of rare causal variants) by the empirical mean (and its corresponding standard error) of the numbers of true discoveries over the 500 simulated datasets divided by the number of causal variants contained in each subgroup.

FDR control

To assess the FDR control of each procedure, we estimated the FDR by the empirical mean (and its corresponding standard error) of the observed false discovery proportion over the 500 simulated datasets.

Overall power

Figure 1 shows the overall power for scenario 1 (reference scenario) with independent markers. Results for scenarios 2 and 3 and for correlated markers are available in supplementary materials. As expected, for all procedures, the overall power tends to decrease with the total number of tested hypotheses m. It also decreases with m₁, since the global effect is distributed among a larger number of SNPs, making the individual effect of each causal variant more difficult to identify. In addition, for correlated datasets, the power tends to increase with ρ in all configurations.

In scenario 1 with independent markers, for small and intermediate values of m₁ (m₁ ≤ 25 and m₁ ≤ 50 for quantitative and binary phenotypes, respectively), wBHa and wBH tend to be the most powerful procedures, although for the smallest m₁ values in the quantitative case (m₁ ≤ 15), BH, qvalue and swfdr are slightly more powerful (Fig. 1). For larger values of m₁, CAMT is the most powerful procedure for quantitative phenotypes (for m₁ ≥ 50) but the least powerful for binary traits (for m₁ ≥ 100), the most powerful being IHW. Conversely, FDRreg has good overall power for binary phenotypes but is the least powerful procedure for quantitative phenotypes. Note that while wBHa is not always the most powerful procedure, it has quite good overall power in all configurations in comparison to the other procedures.

In scenarios 2 and 3 with independent markers, for small and intermediate values of m₁ (m₁ ≤ 50), wBHa, wBH, BH, qvalue and swfdr tend to be the most powerful procedures with similar results (except for binary phenotypes in scenario 3 where BH, qvalue and swfdr are more powerful for m₁ ≤ 50) (Figs. S1 and S2). For larger m₁ values (m₁ ≥ 100), IHW tends to be the most powerful procedure. As in scenario 1, CAMT performs well with quantitative phenotypes but is the least powerful procedure for binary phenotypes.

For correlated markers (Figs. S3, S4 and S5), wBHa is among the most powerful procedures when the value of m₁ is intermediate (in scenario 1) or small (in scenarios 2 and 3). Note that in scenario 1 with binary traits, wBHa tends to be the most powerful procedure in all configurations (Fig. S3).

The results obtained with simulations based on a real dataset are similar to those obtained with fully simulated data (Fig. 2 and Figs. S6 and S7). In all scenarios, for intermediate and small values of m₁ (m₁ ≤ 25), wBHa and wBH tend to be the most powerful procedures while CAMT and IHW are more powerful for large values of m₁ (m₁ ≥ 100).

Power in subgroups

Figure 3 shows the power of the different procedures to detect associations in the subgroup of rare variants for scenario 1 (reference scenario) with independent markers. Results for scenarios 2 and 3 and for correlated markers are available in supplementary materials. As for the overall power, the power in the subgroup of rare variants tends to decrease with the total number of tested hypotheses m and with the number of causal SNPs m₁ for each configuration. In addition, for correlated markers, the power tends to increase with the ρ value in all configurations.

In scenario 1 with independent markers (Fig. 3), the wBH procedure, which is designed particularly for this context, is the most powerful procedure in almost all settings. However, our procedure wBHa shows quite large power compared to the other procedures. For large values of m₁ (m₁ ≥ 50), CAMT tends to be the most powerful procedure for quantitative phenotypes, but it is the least powerful one for binary phenotypes. Conversely, FDRreg performs well for binary phenotypes but it is the least powerful for quantitative phenotypes.

Interestingly, in intermediate scenarios (scenarios 2 and 3), our procedure and wBH tend to be the most powerful in the subgroup of rare variants for all configurations (Figs. S8 and S9). However, as expected, the powers of all procedures for detecting associations of rare variants are low in scenario 2, where the smallest effects are attributed to rare variants.

In the correlation case, for all scenarios we obtained results similar to those obtained in the independent case (Figs. S10, S11 and S12). Thus, wBHa and wBH tend to be the most powerful procedures in the subgroup of rare variants for all settings, except for m₁ ≥ 50 for quantitative traits in scenario 1.

When considering simulations based on a real dataset (Fig. 4 and Figs. S13 and S14), wBHa is among the most powerful procedures in all settings, the most powerful being wBH.

FDR control

Figure 5 displays the estimated FDR for all procedures in scenario 1 with independent markers. These results indicate a good control of the FDR for all procedures in all settings (except for FDRreg with binary phenotypes). Indeed, for all procedures except FDRreg, the estimated FDR is lower than 0.05 or slightly larger than this threshold. This can be explained by the fact that in our simulated settings, the numbers of rejections tend to be small, leading to a quite large variability of the false discovery proportion. Thus, even with the BH procedure (for which the FDR control has been theoretically proven for independent tests), the estimated FDR is slightly larger than the threshold 0.05 in some settings. Similar results were obtained for scenarios 2 and 3 (Figs. S15 and S16).

In the case with correlations between variants (Fig. 6), the estimated FDR increases with ρ in all configurations. Similar results were obtained in scenarios 2 and 3 (Figs. S17 and S18). In addition, the estimated FDR obtained with simulations based on a real dataset (Figs. S19, S20 and S21) tend to be large for all procedures. These results illustrate the difficulty to define and to control the FDR when tested hypotheses are correlated. Nevertheless, as the estimated FDR are similar from one procedure to another, the power comparisons remain relevant.

Larger numbers of tested hypotheses

When considering larger values of m and m₁ for quantitative phenotypes, wBHa belongs to the three most powerful procedures in scenario 1 while maintaining good overall power in scenarios 2 and 3 (Fig. S22). wBHa is also one of the three most powerful procedures in the subgroup of rare variants in scenario 1, while in scenarios 2 and 3, CAMT, IHW and swfdr tend to be slightly more powerful, particularly for large values of m₁ (Fig. S23). Figure S24 indicates a quite good control of the FDR for all procedures.

Other covariate

When MAF is replaced by 1/MAF to prioritize common variants, wBHa is one of the most powerful procedures for detecting common variants in all scenarios (Figs. S25, S26 and S27) while maintaining an overall power similar to that of the unweighted procedure BH (Figs. S28, S29 and S30). The FDR is controlled by all procedures except FDRreg for binary phenotypes (Figs. S31, S32 and S33). Note that the use of log(1/MAF) in CAMT makes it possible to avoid computational problems but leads to very similar results in the cases where CAMT worked with 1/MAF (data not shown).

When using an uninformative covariate, the different procedures remain valid since the FDR is controlled at the desired level (Figs. S34, S35 and S36). All weighted procedures tend to have a lower overall power than the unweighted BH procedure, although wBHa is the one with the smallest loss (Figs. S37, S38 and S39).