GWAR: robust analysis and meta-analysis of genome-wide association studies

Abstract

Motivation

In the context of genome-wide association studies (GWAS), there is a variety of statistical techniques in order to conduct the analysis, but, in most cases, the underlying genetic model is usually unknown. Under these circumstances, the classical Cochran-Armitage trend test (CATT) is suboptimal. Robust procedures that maximize the power and preserve the nominal type I error rate are preferable. Moreover, performing a meta-analysis using robust procedures is of great interest and has never been addressed in the past. The primary goal of this work is to implement several robust methods for analysis and meta-analysis in the statistical package Stata and subsequently to make the software available to the scientific community.

Results

The CATT under a recessive, additive and dominant model of inheritance as well as robust methods based on the Maximum Efficiency Robust Test statistic, the MAX statistic and the MIN2 were implemented in Stata. Concerning MAX and MIN2, we calculated their asymptotic null distributions relying on numerical integration resulting in a great gain in computational time without losing accuracy. All the aforementioned approaches were employed in a fixed or a random effects meta-analysis setting using summary data with weights equal to the reciprocal of the combined cases and controls. Overall, this is the first complete effort to implement procedures for analysis and meta-analysis in GWAS using Stata.

Availability and Implementation

A Stata program and a web-server are freely available for academic users at http://www.compgen.org/tools/GWAR

Supplementary information

Supplementary data are available at Bioinformatics online.

1 Introduction

In genome-wide association studies (GWAS), millions of single nucleotide polymorphisms (SNP) are read in a single genome (Manolio, 2010). Statistical issues in GWAS analysis include quality control, adjustment for population stratification, selection of the optimal statistical test for the analysis and tackling the problem of multiple testing. Some excellent overviews of these issues can be found in (Teo, 2008; Ziegler et al., 2008). Although GWAS have identified many variants associated with complex diseases, these currently explain little of the risk variability for most diseases. Common variants have weak effects and detection of signals requires large sample sizes (Chapman et al., 2011). As single GWAS are often underpowered, meta-analysis, being the statistical procedure for the collection and combination of the available information (Normand, 1999; Petiti, 1994; Trikalinos et al., 2008), increases power and decreases false positive findings. Therefore, most genetic risk variants discovered in the past few years have originated from large-scale meta-analyses of GWAS (Zeggini and Ioannidis, 2009).

In GWAS, the analysis is conducted using a variety of statistical techniques (Balding, 2006; Sasieni, 1997). However, when the genetic model (recessive, additive, dominant and so on) is unknown, robust tests, which retain high power across all plausible models, are preferable (Bagos, 2013; Freidlin et al., 2002). Some of the robust methods are based on the theory of efficiency robust procedures (Freidlin et al., 1999; Gastwirth, 1985), namely the Maximum Efficiency Robust Test (MERT) and the maximum (MAX) test (Freidlin et al., 2002). However, the use of the MAX statistic introduces an extra complexity in the analysis, as under the null hypothesis of no association, it does not follow the standard normal distribution, resulting in computationally intensive resampling-based procedures in order to estimate its P-value. Later, an asymptotic distribution was derived by simplifying the computation of the respective P-values and efficient, faster algorithms were developed (Zang et al., 2010). MAX is uniformly more powerful than Pearson’s chi-square statistic (χ₂²) when the underlying genetic model is recessive, additive, multiplicative and dominant (Joo et al., 2009). In the MIN2 approach, a combination of test statistics is used heuristically and the asymptotic null distribution is derived (Joo et al., 2009). It was shown that MAX has greater efficient robustness than χ₂² and MIN2 when the genetic models are restricted to recessive, additive and dominant (Joo et al., 2009).

Several methods and software tools have been developed for meta-analysis of GWAS [e.g. METAL (Willer et al., 2010) and GWAMA (Magi and Morris, 2010)]. These methods have been employed successfully, with hundreds of GWAS meta-analyses being published (Panagiotou et al., 2013). However, all these tools perform standard meta-analysis assuming the genetic model beforehand (usually, the additive one). A comprehensive comparison of tools is reviewed elsewhere (Evangelou and Ioannidis, 2013). Some of the robust methods allow for a stratified analysis, which can be formulated as a fixed-effect meta-analysis (Hothorn and Hothorn, 2009; Pan et al., 2011; Zang and Fung, 2010, 2011; Zheng and Tian, 2006), whereas other approaches are direct extensions of the traditional methods of meta-analysis under random-effects models (Pereira et al., 2011; Zintzaras, 2010). Another approach is to perform a separate analysis for each study, using the robust method of choice (i.e. MAX, MERT or MIN2) and then combine the individual results by pooling the Z-values or the effect sizes in a fixed or random effects method, as described by Zhou et al. (2011). This approach, which is closer to the norm followed for meta-analysis of GWAS (de Bakker et al., 2008), requires only summary estimates and was adopted in this work.

The objective of this work is to present efficient implementations of the Cochran-Armitage trend test (CATT) under a recessive, additive and dominant model of inheritance and, most importantly, of the robust methods MERT (Freidlin et al., 2002), MAX (Zang et al., 2010) and MIN2 (Joo et al., 2009) in Stata and to offer an easy-to-use, fast and reliable software. All the aforementioned approaches were also employed in a fixed or a random effects meta-analysis setting for summary data.

2 Materials and methods

2.1 Notation

Consider a case-control association study with a biallelic locus (i.e. a SNP) with alleles A and B having frequencies P_A and P_B, respectively. The three possible genotypes are denoted by g₀ = AA, g₁ = AB and g₂ = BB with respective frequencies π_j = P(g_j) for j = 0,1,2. We further denote the disease penetrances by f_j = P(case|g_j) and the disease prevalence by k = P(case) = Σf_jπ_j. In a case-control study, r cases and s controls are sampled independently and the observed genotype counts (r₀, r₁, r₂) and (s₀, s₁, s₂) follow multinomial distributions, mul(r; p₀, p₁, p₂) and mul(s; q₀, q₁, q₂), respectively (Table 1), with p_j = P(g_j|case) = f_jπ_j/k and q_j = P(g_j|control) = (1−f_j)π_j/(1−k) for j = 0,1,2. The null hypothesis is equivalent to H₀: f₀ = f₁ =  f₂ = k. Assuming that B is the risk allele, we may denote the genotype relative risks by γ_j = f_j/f₀ for i = 1,2 and f₀ > 0. Then, the null hypothesis can be reformulated as H₀: γ₂ = γ₁ = 1 against H₁: γ₂ ≥ γ₁ ≥ 1 with at least one inequality. Standard models for disease penetrance that imply a specific relationship between genotype and disease include common dominant and common recessive models, multiplicative or additive models. Thus, the underlying genetic model of inheritance can be dominant, recessive, multiplicative or additive, if γ₂ = γ₁ > 1, γ₂ > γ₁ = 1, γ₂ = 2γ₁−1 > 1 and γ₂ = γ₁² > 1, respectively. The multiplicative and additive models are asymptotically equivalent and in many cases are indistinguishable (hence the general term co-dominant may be used sometimes) (Bagos, 2013; Clarke et al., 2011). Another less frequently encountered situation occurs when γ₁ > γ₂ = 1, in which case the genetic model is termed overdominant.

Table 1

The 2 × 3 contingency table with the distribution of cases and controls in a traditional GWAS concerning a single biallelic locus

	AA (g0)	AB (g1)	BB (g2)	Total
Cases	r0	r1	r2	r
Controls	s0	s1	s2	s
Total	n0	n1	n2	n

2.2 Analysis of GWAS

There is a variety of methods for analyzing data from GWAS (Balding, 2006; Langefeld and Fingerlin, 2007). Perhaps the most common analysis of genotypes in case-control studies under a genetic model-free (GMF) perspective is to examine the association between the rows and columns of the 2 × 3 contingency table (Table 1). This can be performed using the traditional Pearson’s chi-square statistic (χ₂²), which is distributed following a chi-square distribution with two degrees of freedom (Balding, 2006; Langefeld and Fingerlin, 2007). Another option would be to calculate the two odds-ratios that correspond to the comparison of the genotypes carrying the risk allele (i.e. g₁ and g₂) against g₀. Such tests can be performed in a logistic regression framework using the case/control status as the dependent variable (Balding, 2006). Alternatively, one may compute the odds-ratios using the summary counts of the contingency table (Chen and Chatterjee, 2007). Asymptotically, these methods are equivalent to χ₂², but the major disadvantage of all these GMF methods lies in the lack of power.

At the same time, methods that are based on prior assumption on the genetic model also exist. These methods involve some sort of re-arrangement of the data or a differential scoring of the genotypes (Sasieni, 1997; Zheng et al., 2003). For instance, if we re-arrange the data of Table 1 by grouping genotypes g₁ and g₂ and perform a χ₁², we explicitly assume a dominant model of inheritance; if we group g₀ and g₁ and compare them to g₂, we are testing the recessive model. An allele-based comparison implies an additive genetic model under the assumption of HWE. Nevertheless, the most commonly used test is CATT (Armitage, 1955; Cochran, 1954), which is a score test for linear trend of the proportions in the contingency table. Assuming a general scoring scheme x_i = (x₀, x₁, x₂) for the genotypes (g₀, g₁, g₂), a Z-test is formulated as:

$$\begin{array}{l}{Z}_{CATT\left(x\right)}=\frac{U_x}{\sqrt{{\mathrm{var}}_{H_0}\left(U_x\right)}}\\ =\frac{\sqrt{n}{\sum }_{i=0}^2{x}_i\left(s{r}_i-r{s}_i\right)}{\sqrt{rsn\left[n{\sum }_{i=0}^2{x}_i^2{n}_i-{\left({\sum }_{i=0}^2{x}_i{n}_i\right)}^2\right]}}\sim N(0,1)\end{array}$$

(1)

where U_x is the score statistic and ${\mathrm{var}}_{Ho}(U_x)$ its variance evaluated under H₀. The scores are usually set to x_i = (0, ½, 1), in which case the statistic tests the additive model of inheritance (Z_CATT(½)). CATT is invariant under linear transformations and, additionally, it does not rely on the assumption of HWE. It should be noted that the corresponding CATT is the most powerful statistic under the dominant, the additive and the recessive model, a fact that explains its widespread use in GWAS (Zheng et al., 2003). The CATT(½) is routinely used in real life applications, since it is more powerful when compared against the χ₂² and does not assume HWE.

2.3 Robust analysis of GWAS

The most prominent methods for robust analysis are based on the theory of efficiency robust procedures (Freidlin et al., 1999; Gastwirth, 1985), namely the MERT and the MAX test (Freidlin et al., 2002). For a family of optimal tests, MERT is a linear combination of the two optimal tests for the extreme members of the family, i.e., the tests with the minimum correlation (ρ). Two conditions need to be satisfied; (i) ρ ≥ 0, and (ii) ρ + 1 must be less than the sum of the other two correlations. Henceforth, we focus on the optimal tests for recessive, additive and dominant models. Freidlin et al. (2002) showed that conditions (i) and (ii) hold if recessive and dominant relationships are modeled. Thus, MERT is expressed as a linear combination of the optimal CATTs for the recessive and dominant models and is asymptotically normally distributed:

$$Z_{MERT}=\frac{Z_{CATT\left(0\right)}+Z_{CATT\left(1\right)}}{\sqrt{2\left(1+{\rho }_{CATT\left(0,1\right)}\right)}}\sim N\left(0,1\right)$$

(2)

In Eq. (2), Z_CATT(0) and Z_CATT(1) denote the Z-values under a recessive and dominant model respectively and ρ_CATT(0,1) denotes the correlation for (Z_CATT(0), Z_CATT(1)) (Freidlin et al., 2002). Expressions for the correlations among the three trend tests were correctly presented in (Zang et al., 2010). The MAX test is based on the simple idea to test all three possible models and choose the one with the highest score (note that we use the absolute value for completeness, in case A is the risk allele):

$$Z_{MAX}=\max (\left|Z_{CATT(0)}\right|,\left|Z_{CATT(1/2)}\right|,\left|Z_{CATT(1)}\right|)$$

(3)

However, for calculating the P-value (i.e. P(|Z_MAX|>t)), we need to take into account the correlations in order to avoid inflation of type I error rate. Freidlin and coworkers estimated the correlations of the three test statistics under H₀ (Freidlin et al., 2002), but, even in that case, simulations or numerical integration are needed; this makes the method computationally intensive. MAX is more powerful, but MERT is easier to apply since the approximate P-value can be obtained from a standard normal distribution (Freidlin et al., 2002). Traditionally, Monte Carlo simulations were used for P-value calculations of MAX, but several less computationally intensive approximations (rhombus formula) have been described in Li et al. (2008). Recently, Zang et al. found that Z_CATT(0), Z_CATT(½) and Z_CATT(1) are linearly dependent, a result that allowed them to develop faster algorithms for calculating the statistical significance of MAX (Zang et al., 2010). Thus, the P-value for the MAX statistic is given by:

$$\begin{array}{l}P\left(Z_{MAX}<t\right)=2{\int }_0^{t\left(1-{\omega }_1\right)/{\omega }_0}\Phi \left(\frac{t-{\rho }_{CATT\left(0,1\right)}{z}_0}{\sqrt{1-{\rho }_{CAT{T}_{\left(0,1\right)}^2}}}\right)\varphi \left(z_0\right)d{z}_0\\ +2{\int }_{t\left(1-{\omega }_1\right)/{\omega }_0}^t\Phi \left(\frac{\left(t-{\omega }_0{z}_0\right)/{\omega }_1-{\rho }_{CATT\left(0,1\right)}{z}_0}{\sqrt{1-{\rho }_{CAT{T}_{\left(0,1\right)}^2}}}\right)\varphi \left(z_0\right)d{z}_0\\ -2{\int }_0^t\Phi \left(\frac{-t-{\rho }_{CATT\left(0,1\right)}{z}_0}{\sqrt{1-{\rho }_{CAT{T}_{\left(0,1\right)}^2}}}\right)\varphi \left(z_0\right)d{z}_0\end{array}$$

(4)

where:

$${\omega }_0=\frac{{\rho }_{CATT\left(0,1/2\right)}-{\rho }_{CATT\left(0,1\right)}{\rho }_{CATT\left(1/2,1\right)}}{1-{\rho }_{CAT{T}_{\left(0,1\right)}^2}}$$

(5)

$${\omega }_1=\frac{{\rho }_{CATT\left(1/2,1\right)}-{\rho }_{CATT\left(0,1\right)}{\rho }_{CATT\left(0,1/2\right)}}{1-{\rho }_{CAT{T}_{\left(0,1\right)}^2}}$$

(6)

In Equation (5–6), ρ_CATT(0,1/2) denotes the correlation (Z_CATT(0), Z_CATT(1/2)) and ρ_CATT(1/2,1) the correlation (Z_CATT(1/2), Z_CATT(1)). MIN2 is an interesting robust approach that was adopted by investigators of the Wellcome Trust Case-Control Consortium (WTCCC, 2007). They applied the χ₂² along with the CATT(½) and, subsequently, chose the minimum of the P-values:

(7)

The use of MIN2 is justified by simulations showing that χ₂² has 5% less power compared with MAX and outperforms MERT, except when the additive model holds (Zheng et al., 2006). However, MIN2 is not a proper P-value since the statistics are correlated and multiple tests are performed. Later, Joo et al. (2009) derived the joint distribution needed in order to calculate a proper P-value:

$$\begin{array}{l}P\left(Z_{CATT\left(1/2\right)}^2<t_1,T_{{\chi }_2^2}^2<t_2\right)=\\ \{\begin{array}{l}1-\frac{1}{2}{e}^{-\frac{t_{{}_1}}{2}}-\frac{1}{2}{e}^{-\frac{t_2}{2}}+\frac{1}{2\pi }{\displaystyle {\int }_{t_1}^{t_2}{e}^{-\frac{u}{2}}}\arcsin \left(\frac{2t_1}{u}-1\right)du,t_1<t_2\\ 1-e^{-\frac{t_2}{2}},t_1>t_2\end{array}\end{array}$$

(8)

Unlike MAX, MIN2 is independent of the allele frequency. Additionally, the P-value of MIN2 is strictly increasing in MIN2. Thus, MIN2 can be directly used to rank all SNPs as WTCCC (2007) did. Simulations have shown that MIN2 is slightly less powerful (1–4%) than MAX when the recessive model holds; it outperforms MAX under the additive model, whereas, for the dominant model, both tests perform similarly. We also observe that MAX needs three integrals whereas MIN2 only one and thus it is slower. The actual complexity depends on the number of integration points used (in our implementation we chose 100).

2.4 Meta-analysis of GWAS

Meta-analysis is a method that aims at increasing the power to detect weak genotype effects, estimating the effect size with greater precision and investigating between-studies variability. Meta-analysis of GWAS using summary statistics, in contrast to direct analysis of pooled individual-level data, diminishes the limitations that are imposed by restrictions on sharing individual-level data. Empirical findings (Janssens et al., 2009) suggest that using time-consuming individual data methods for meta-analysis of GWAS does not necessarily lead to increase efficiency (Lin and Zeng, 2010a,b) and thus, summary methods are routinely used for meta-analyses of GWAS (de Bakker et al., 2008; Magi and Morris, 2010; Willer et al., 2010).

Some of the statistics discussed in the previous sections, such as the CATT, are formed as a score test:

$$Z_{score\left(i\right)}=\frac{U_{{\beta }_i}}{\sqrt{{\mathrm{var}}_{H_0}\left(U_{{\beta }_i}\right)}}=\frac{U_{{\beta }_i}}{\sqrt{\widehat{I}\left({\beta }_i\right)}}$$

(9)

where U_βi is the score statistic from study i, and $\widehat{I}({\beta }_i)$ the estimate of the information matrix evaluated at H₀ (i.e. $\widehat{{\beta }_i}$=0). Other tests (i.e. those based on the log-odds ratios) are formulated as a Wald statistic:

$$Z_{Wald\left(i\right)}=\frac{{\widehat{\beta }}_i}{\sqrt{\mathrm{var}\left({\widehat{\beta }}_i\right)}}=\frac{{\widehat{\beta }}_i}{\sqrt{\widehat{I}{\left({\beta }_i\right)}^{-1}}}={\widehat{\beta }}_i\sqrt{\widehat{I}\left({\beta }_i\right)}$$

(10)

Here, ${\widehat{\beta }}_i$ is the maximum likelihood estimate and its variance $\mathrm{var}({\widehat{\beta }}_i)$ is calculated at ${\widehat{\beta }}_i$⁠. If both tests evaluate the same null hypothesis and if ${\widehat{\beta }}_i$ is small in absolute value, the score and Wald statistics are asymptotically equivalent under H₀ (i.e. ${\widehat{\beta }}_i$ = 0).

Even though optimal methods for combining the results of independent studies that used a Wald statistic are available for years and described in standard textbooks (Normand, 1999; Petiti, 1994), optimal weights for combining the results of score statistics such as the CATT had not been presented until recently (Zhou et al., 2011). Consider a meta-analysis of k studies under a fixed-effect model in which the true effect is common across the k studies (i.e. β_i ∼ N(β, s_i²)). Then, the pooled Z will be given by the weighted average:

$$\overline{Z}=\frac{{\displaystyle {\sum }_{i=1}^k\sqrt{w_i}{Z}_i}}{\sqrt{{\displaystyle {\sum }_{i=1}^k{w}_i}}}$$

(11)

Zhou and coworkers have shown that, in any case, the optimal weights will be given by $w_i=\widehat{I}({\beta }_i)$ (Zhou et al., 2011). It should be noted that when we are pooling Z-values, irrespective of the method by which they have arisen, $\widehat{I}({\beta }_i)$ is usually not available to the meta-analyst. An approach commonly followed in the literature is to weigh the studies proportionally to the study-specific sample sizes (de Bakker et al., 2008; Magi and Morris, 2010; Willer et al., 2010). However, Zhou et al. showed that even though this scheme performs well for most practical applications, it is suboptimal and proposed that a good choice for the weights would be $w_i=r_i{s}_i/n_i={(1/r_i+1/s_i)}^{-1}$⁠. Equation (11) is a general formula and in the case of pooling a score statistic, it reduces to:

$${\overline{Z}}_{score}=\frac{{\sum }_{i=1}^k{U}_{{\beta }_i}}{\sqrt{{\sum }_{i=1}^k\widehat{I}\left({\beta }_i\right)}}$$

(12)

which is a well-known expression for the stratified score statistic obtained from the sum of k log-likelihoods (i.e. a Mantel-Haenszel type statistic). In case we pool a Wald statistic, Equation (11) is equivalent to obtaining a pooled estimate of β:

$$\widehat{\beta }=\frac{{\displaystyle {\sum }_{i=1}^k{w}_i{\beta }_i}}{{\displaystyle {\sum }_{i=1}^k{w}_i}}$$

(13)

This estimate is the well-known inverse-variance estimate commonly used in meta-analysis, which utilizes the estimated standard errors of the beta coefficients as weights. Asymptotically, weighting by the sample size (or the method proposed by Zhou and coworkers) and the standard errors are equivalent when the trait distribution is identical across samples (such that standard errors are a predictable function of sample size).

In the presence of between-studies heterogeneity, a random-effects model for pooling the estimates is preferable, assuming β_i ∼ N(β, s_i² + τ²) with τ² being the between-studies variance estimate. Usually, τ² is estimated from the data using the method of moments, even though there are several alternatives. The problem, however, with robust statistics is that only a P-value or a Z-score is reported, and not a Wald statistics along with its variance. Thus, we relied on an approximation proposed by Zhou and coworkers (Zhou et al., 2011). They used the same method for modeling the components of the score test and showed that, when β_i is small, the asymptotic mean of U_βi can be approximated by ${\beta }_i\widehat{I}({\beta }_i)$⁠. Thus, after calculating the appropriate weights equal to $w_i=\widehat{I}({\beta }_i)={(1/r_i+1/s_i)}^{-1}$⁠, we use the Z-score to calculate an effect size, taking the direction of the association into account:

$${\widehat{\beta }}_i=\frac{U_{{\beta }_i}}{\widehat{I}\left({\beta }_i\right)}=\frac{Z_i}{\sqrt{\widehat{I}\left({\beta }_i\right)}}=\frac{Z_i}{\sqrt{{\left(1/r_i+1/s_i\right)}^{-1}}}$$

(14)

We can then perform a standard random-effects meta-analysis; we only need to keep in mind that the actual effect sizes calculated by Equation (14), as well as by the meta-analysis, are not valid. Consequently, we have to rely only on the overall Z-score and the overall P-value, which are valid under any circumstances. This method can be regarded as an extension of the Stouffer’s method (Equation 11)), which was one of the first methods that combined P-values. Thus, in principle other methods that combine P-values could be used, but these have serious disadvantages: the direction of the association is not taken into consideration and therefore all P-values have to be one-sided, and, most importantly, they do not account for between-studies heterogeneity.

3 Results

3.1 Implementation

GWAR is implemented in Stata (StataCorp). For the MAX and MIN2 statistics, we implement the methods described in Equations (4) and (8), respectively, that require integration. For this purpose, we use numerical integration with Gauss-Legendre quadrature implemented in the integrate command that uses Mata, and this results in great gain in computing time. For the meta-analysis part, we use the metan command in Stata. Due to integrate, GWAR requires Stata v12 or above. GWAR is implemented as an ado file as well as a web-server available at www.compgen.org/tools/GWAR.

3.2 Usage

There are two programs: gwar and min2gwar. The syntax for gwar is: gwar varlist, Method(string) [SNP(varlist min = 1 max = 1) EFfect(string)]. varlist contains the variables corresponding to the genotypes for cases and controls. The first three variables are the genotypes of controls (e.g. aa0, ab0, bb0), and the remaining ones the genotypes of cases (e.g. aa1, ab1, bb1). Allele b is assumed to be the risk variant. In the single-study mode, the user needs to specify only the genotypes and the method of analysis (MAX, MIN2, MERT as well as the three versions of CATT). In the meta-analysis mode, the user also needs to specify the variable that indicates the SNPs and the requested method of meta-analysis (fixed or random effects). Meta-analysis is performed, using Stata’s metan command, by grouping the observations with the same SNP, so it is the user’s responsibility to ensure that data are in the appropriate format when the program is called. For fast grouping of large datasets the program uses the ftools library, which the user has to install.min2gwar is the version that uses the P-values corresponding to the Pearson′s chi-square test and the CATT(1/2) directly. In the single-study mode, its use is simple (min2gwar p1 p2). In the meta-analysis mode, the user needs also to provide three additional variables, the total number of cases and controls and a variable (−1, +1) for the direction of the association. Similarly, the meta-analysis is performed by grouping the observations with the same SNP, so it is the user's responsibility to ensure that data are in the appropriate format when the program is called.

The web-server provides the functionality for uploading a file, or it can be used as a simple online calculator in which the user can fill the genotype counts and receive the results. For large datasets, we allow uploading compressed files in most popular formats, and we also allow the user to receive the results by email.

3.3 Performance

GWAR is fast, reliable and possesses features currently not available in other similar packages. The benchmark was performed on an Intel Xeon CPU E5-2660v2, 2.20GHz server with 32 GB of RAM using STATA/SE. A GWAS with 100 000 SNPs can be analyzed within 1-2 hours, even using the computationally intensive MAX or MIN2 methods (the analysis using CATT or MERT is performed within seconds), whereas, for a meta-analysis of five studies and 100 000 SNPs, the results are produced depending on the method, from half an hour to four hours (Supplementary Table S1). The accuracy can be seen with a comparison against results known in the literature. We need to mention that, due to the approximate formula employed, GWAR can calculate P-values up to 10⁻¹⁸ for MAX and 10⁻⁵² for MIN2, values that are nevertheless outside of the range of any reasonable simulation method (for the latter one would need more than 10⁵² repetitions, a number which is untenable).

In order to compare the robust methods (MERT, MAX and MIN2), we performed a simulation study. We simulated 500 cases and 500 controls. We considered different values for the risk allele frequency (RAF) ranging from 0.1 to 0.5, the odds-ratio for the g₂ versus g₀ comparison (1–2), and the genetic model of inheritance (dominant, additive and recessive). The results indicate that, for each mode of inheritance, the respective CATT is, as one would expect, more powerful compared with every other test statistic (see Supplementary Table S1, Supplementary Fig. S1). For the recessive model, MAX is more powerful than MERT and MIN2, with differences ranging from 0.2 to 13.2%. Moreover, MAX is almost uniformly more powerful than CATT(1/2) which is routinely used in GWAS with a difference in power that exceeds 10% for odds-ratios >1.5. The latter also holds for MERT which can be used as an alternative to reduce complexity that arises from the computation of MAX and MIN2. Under the additive model, MIN2 outperforms MAX, but the difference in power is <3% for odds-ratios <1.5. For the dominant model, MAX and MIN2 have similar performance. MERT is the least powerful, with a difference in power up to 17.4% for RAF = 0.1 when compared against MAX.

3.4 Application to real studies

For illustration purposes, we re-analyzed the 22 SNPs with an estimated P-value [by the rhombus formula for approximating the P-value of MAX (Li et al., 2008)] below 10⁻⁵ from a GWAS of Coronary Artery Disease (CAD) (WTCCC, 2007). The results from the MAX method of analysis in GWAR have a remarkable agreement with those obtained using the rhombus formula (Fig. 1). It should also be noted that both approaches can be used to estimate extremely low P-values (say <10⁻⁷) where resampling procedures are not easily applicable (or they are too demanding).

Fig. 1

The comparison of P-values obtained from the MAX method of analysis in GWAR versus those obtained using the rhombus formula for 22 chosen SNPs from the GWAS of CAD

We also analyzed the data from a GWAS that was designed to identify susceptibility loci for type-2 diabetes (Sladek et al., 2007). We calculated the P-values from the MAX method of analysis in GWAR and compared them to those obtained by the authors using a permutation approach as well as using a proposed asymptotic distribution of the MAX statistic (Gonzalez et al., 2008). We observe an almost perfect correlation among the results with consistently, but only slightly higher P-values obtained by GWAR (Fig. 2).

Fig. 2

The comparison of P-values obtained from the MAX method of analysis in GWAR versus those obtained using a permutation approach and the asymptotic distribution of the MAX statistic (Gonzalez et al., 2008) for 67 SNPs from the GWAS of type-2 diabetes. Four extra SNPs (rs7903146, rs12255372, rs10885409 and rs7904519) were also reported, all with estimated P-values: <3.30E-10 for the permutation approach, 6.36E-19, 9.33E-15, 2.47E-11 and 5.60E-11 for the asymptotic procedure and 3.77E-15, 1.89E-14, 4.21E-11 and 6.44E-11 for GWAR, respectively

We finally analyzed 6 SNPs from a hypertension GWAS (WTCCC, 2007). We calculated the P-values from the MAX and MIN2 method of analysis in GWAR and compared them to those reported by (Joo et al., 2009). The results are identical for MIN2 with some small discrepancies for the MAX test (Table 2).

Table 2

Application to the SNPs associated with hypertension in WTCCC (2007) with MIN2 and MAX robust tests

	MIN2 reported in (Joo et al., 2009)	MIN2 in gwar	MAX reported in (Joo et al., 2009)	MAX in gwar
rs2820037	1.31E-06	1.31E-06	3.16E-06	3.23E-06
rs6997709	1.34E-05	1.34E-05	2.06E-05	2.07E-05
rs7961152	1.25E-05	1.25E-05	1.80E-05	2.01E-05
rs11110912	1.56E-05	1.56E-05	8.16E-06	8.15E-06
rs1937506	1.56E-05	1.56E-05	2.46E-05	2.43E-05
rs2398162	9.64E-06	9.64E-06	2.50E-06	2.42E-06

4 Discussion

We presented GWAR, a fast and reliable tool for analysis and meta-analysis of GWAS. We implemented both standard CATTs as well as the most widely used robust methods, namely MERT, MAX and MIN2. For MAX and MIN2, which are also the most powerful among the robust methods, we implemented accurate and fast routines that rely on numerical integration in order to calculate the P-value accurately. GWAR is implemented as an ado file for Stata 12 or higher, as well as a web-server available at www.compgen.org/tools/GWAR.

GWAR offers functionalities not available in other similar tools for the analysis of GWAS (coin (Hothorn and Hothorn, 2009), SNPassoc (Gonzalez et al., 2007), Rassoc (Zang et al., 2010), RobustSNP (So and Sham, 2011), GenABEL (Aulchenko et al., 2007) and Plink (Purcell et al., 2007). First of all, it implements MIN2, which is currently not available in any statistical package. MAX is also available in coin, SNPassoc, Rassoc and RobustSNP, whereas MERT is available only in Rassoc. Even though MAX implementation in Stata is among the second slowest found in these tools (it surpasses only RobustSNP), MERT and the various CATTs are the fastest ones (Supplementary Tables S1 and S3).

GWAR is also unique in terms of meta-analysis since it is the only tool that implements meta-analysis methods using the aforementioned robust statistics. Similar methods are not available in other popular packages for meta-analysis of GWAS, like METAL (Willer et al., 2010), GWAMA (Magi and Morris, 2010), PLINK (Purcell et al., 2007) and MetABEL (Aulchenko et al., 2007). Even for the standard meta-analysis that assumes a model of inheritance beforehand, GWAR is unique since it offers the option for both fixed and random effects meta-analysis. For comparison purposes, we need to mention that only PLINK and GWAMA offer the option for random effects meta-analysis, whereas MetABEL and METAL allow only fixed effect. Moreover, GWAMA, PLINK and METAL, although computational optimal, require precalculated log-odds ratios from the individual studies assuming a specific model of inheritance beforehand. PLINK and METAL offer also an option for Z-score-based P-values meta-analysis, weighted by the sample size. Finally, MetABEL allows for meta-analysis of CATT(1/2) using the qtscore function.

The web-server that we implemented is also unique of its kind since it allows researchers to upload files and receive the results without having to run the software locally. The tool requires only summary genotype counts, and thus security is not a main concern (even though we guarantee the privacy of the results, and in any case the user may not provide the list with the actual SNP ids). Nevertheless, there is theoretical evidence that summary statistics provide the same efficiency of the estimates (Lin and Zeng, 2010a,b) and meta-analysis of GWAS is usually performed using summary estimates from the included studies (Begum et al., 2012; de Bakker et al., 2008; Magi and Morris, 2010; Willer et al., 2010). However, MIN2 statistic can also be used even without the genotype counts. The user may only provide the P-values from the Pearson’s chi-square and the standard CATT (for the single study analysis), or additionally the total number of cases and controls and the direction of the association (for the meta-analysis). With this option, the analyst may need to perform some calculations, but he/she avoids additional complications arising from confidentiality agreements, especially when it comes to large consortia in meta-analysis of GWAS.

Our simulations showed that MAX is more powerful than MIN2 (except when the additive model holds), and MIN2 is always more powerful than MERT. However, the computation time needed is in the reverse order (MERT is faster than MIN2 and MIN2 faster than MAX). All these should be taken into account by the user. For instance, if time is not an issue, MAX should be used. MIN2 is a good compromise (especially if we consider that most associations behave like additive), whereas MERT is fast and could be used for really large datasets, or as an exploratory tool.

Currently, GWAR allows only for analysis of case-control studies. Possible future improvements of GWAR will include the ability to analyse also quantitative traits. To this end however, only the MERT and MAX statistics will be possible to implement, since MIN2 is inherently bound to case-control studies. Another possible extension would be to allow for covariates since it is common that GWAS suffer from several hidden biases that need to be adjusted for. In most GWAS, some type of adjustment is used by allowing for covariates (e.g. gender, age, eigenvectors), but the majority of robust methods are not easily amenable for such modifications. In principle, this could be possible for GWAR by modifying the score statistic and using individual data. This would be a major update of the tool and we plan it for the future. However, one of the most important sources of bias in GWAS, which is population stratification, that can inflate the test statistics, can be estimated and adjusted for, without the modification of the basic method (Clayton et al., 2005). One can generate a quantile–quantile plot of the observed test statistics for disease association versus each of the SNPs included in the analysis. Usually, in such approaches SNPs with extreme deviations from HWE in either controls or cases are excluded, and the same has to be done for the top 10% of the SNPs. The tests can then be ranked from smallest to largest and plotted against their expected values under the global null hypothesis that no true association exists. The slope of the fitted line would produce an estimate of the inflation and suggest adjustments to the significance level. Finally, another possible extension will be to implement the tool in a standalone version (e.g. in Perl programming language), in order to increase the speed and minimize the dependence on external software.

Acknowledgements

The authors also gracefully acknowledge NGP-net COST Action (BM1405) for support and the three reviewers for their valuable comments and suggestions that improved the quality of the article.

Funding

The authors are financially supported by the project ‘Integration of Data from Multiple Sources’ (IntDaMus), which is implemented under the ‘ARISTEIA ΙΙ’ Action of the ‘OPERATIONAL PROGRAMME EDUCATION AND LIFELONG LEARNING’ and is co-funded by the European Social Fund and National Resources.

Conflict of Interest: none declared.