Gene analysis

The gene analysis in MAGMA is based on a multiple linear principal components regression [18] model, using an F-test to compute the gene p-value. This model first projects the SNP matrix for a gene onto its principal components (PC), pruning away PCs with very small eigenvalues, and then uses those PCs as predictors for the phenotype in the linear regression model. This improves power by removing redundant parameters, and guarantees that the model is identifiable in the presence of highly collinear SNPs. By default only 0.1% of the variance in the SNP data matrix is pruned away.

With $X_g^{*}$ the matrix of PCs, Y the phenotype and W an optional matrix of covariates the model can thus be written as $Y={\alpha }_{0g}\overrightarrow{1}+X_g^{*}{\alpha }_g+W{\beta }_g+{\epsilon }_g$, where the parameter vector α _g represents the genetic effect, β _g the effect of the optional covariates, α _0g the intercept and ε _g the vector of residuals. The F-test uses the null-hypothesis H ₀: ${\alpha }_g=\overrightarrow{0}$of no effect of gene g on the phenotype Y, conditional on all covariates.

This choice of gene analysis model is motivated by a balance of statistical and practical concerns. This multiple regression model ensures that LD between SNPs is fully accounted for. It also offers the flexibility to accommodate additional covariates and interaction terms as needed without changing the model. At the same time, since the F-test has a known asymptotic sampling distribution the gene p-values take very little time to compute, making the gene analysis much faster than permutation-based alternatives.

The linear regression model is also applied when Y is a binary phenotype. Although this violates some assumptions of the F-test, comparison of the F-test p-values with p-values based on permutation of the F-statistic shows that the F-test remains accurate (see ‘Supplemental Methods—Implementation Details’). MAGMA therefore uses the asymptotic F-test p-values by default, though it also offers an option to compute permutation-based p-values using an adaptive permutation procedure. In addition, comparison with logistic regression models shows that the results of the linear model are effectively equivalent to that of the more conventional logistic regression model, but without the computational cost.

Gene-set analysis

To perform the gene-set analysis, for each gene g the gene p-value p _g computed with the gene analysis is converted to a Z-value z _g = Φ⁻¹(1 – p _g), where Φ⁻¹ is the probit function. This yields a roughly normally distributed variable Z with elements z _g that reflects the strength of the association each gene has with the phenotype, with higher values corresponding to stronger associations.

Self-contained gene-set analysis tests whether the genes in a gene-set are jointly associated with the phenotype of interest. As such, using this variable Z a very simple intercept-only linear regression model can now be formulated for each gene set s of the form$Z_s={\beta }_0\overrightarrow{1}+{\epsilon }_s$, where Z _s is the subvector of Z corresponding to the genes in s. Evaluating β ₀ = 0 against the alternative β ₀ > 0 yields a self-contained test, since under the self-contained null hypothesis that none of the genes is associated with the phenotype z _g has a standard normal distribution for every gene g.

Competitive gene-set analysis tests whether the genes in a gene-set are more strongly associated with the phenotype of interest than other genes. To test this within the regression framework the model is first expanded to include all genes in the data. A binary indicator variable S _s with elements s _g is then defined, with s _g = 1 for each gene g in gene set s and 0 otherwise. Adding S _s as a predictor of Z yields the model$Z={\beta }_{0s}\overrightarrow{1}+S_s{\beta }_s+\epsilon$. The parameter β _s in this model reflects the difference in association between genes in the gene set and genes outside the gene set, and consequently testing the null hypothesis β _s = 0 against the one-sided alternative β _s > 0 provides a competitive test. Note that this is equivalent to a one-sided two-sample t-test comparing the mean association of gene-set genes with the mean association of genes not in the gene-set. Similarly, the self-contained analysis is equivalent to a one-sided single-sample t-test comparing the mean association of gene-set genes to 0.

It should be clear that in this framework, the gene-set analysis models are a specific instance of a more general gene-level regression model of the form$Z={\beta }_0\overrightarrow{1}+C_1{\beta }_1+C_2{\beta }_2+\dots +\epsilon$. The variables C ₁, C ₂, …, in this generalized gene-set analysis model can reflect any gene property, from the binary indicators used for the competitive gene-set analysis to continuous variables such as gene size and expression levels. Any transformations of, and interactions between, such gene properties can also be added. This generalized gene-set analysis model thus allows for testing of conditional, joint and interaction effects of any combination of gene sets and other gene properties. In practice, the competitive gene-set analysis implemented in MAGMA in fact uses such a generalized model by default, performing a conditional test of β _s corrected for the potentially confounding effects of gene size, gene density and (if applicable, e.g. in meta-analysis) difference in underlying sample size, if such effects are present. This is achieved by adding these variables, as well as the log of these variables, as covariates to the gene-level regression model. The gene density is defined as the ratio of effective gene size to the total number of SNPs in the gene, with the effective gene size in turn defined as the number of principal components that remain after pruning.

One complication that arises in this gene-level regression framework is that the standard linear regression model assumes that the error terms have independent normal distributions, i.e.$\epsilon ~\text{MVN}(\overrightarrow{0},{\sigma }^{2}I)$. However, due to LD, neighbouring genes will generally be correlated, violating this assumption. This issue can be addressed by using Generalized Least Squares approach instead, and assuming that $\epsilon ~\text{MVN}(\overrightarrow{0},{\sigma }^{2}R)$. In MAGMA, the required gene-gene correlation matrix R is approximated by using the correlations between the model sum of squares (SSM) of each pair of genes from the gene analysis multiple regression model, under their joint null hypothesis of no association. These correlations are a function of the correlations between the SNPs in each pair of genes and thus provide a good reflection of the LD, and since they have a convenient closed-form solution they are easy to compute (see also ‘Supplemental Methods—Implementation Details’). Note that for the self-contained analysis, the submatrix R _s corresponding to only the genes in the gene set is used instead of R. In addition, since the self-contained null hypothesis guarantees that all z _g have a standard normal distribution, the error variance σ ² can be set to 1.

Analysis of summary SNP statistics

Since raw genotype data may not always be available for analysis, MAGMA also provides more traditional SNP-wise gene analysis models of the type implemented in PLINK and VEGAS. These SNP-wise models first analyse the individual SNPs in a gene and combine the resulting SNP p-values into a gene test-statistic, and can thus be used even when only the SNP p-values are available. Although evaluation of the gene test-statistic does require an estimate of the LD between SNPs in the gene, estimates based on reference data with similar ancestry as the data the SNP p-values were computed from has been shown to yield accurate results [17,19].

Two types of gene test statistics have been implemented in MAGMA: the mean of the χ ² statistic for the SNPs in a gene, and the top χ ² statistic among the SNPs in a gene. For the mean χ ² statistic, a gene p-value is then obtained by using a known approximation of the sampling distribution [20,21]. For the top χ ² statistic such an approximation is not available, and therefore an adaptive permutation procedure is used to obtain an empirical gene p-value. A random phenotype is first generated for the reference data, drawing from the standard normal distribution. This is then permuted, and for each permutation the top χ ² statistic is computed for every gene. The empirical p-value for a gene is then computed as the proportion of permuted top χ ² statistics for that gene that are higher than its observed top χ ² statistic. The required number of permutations is determined adaptively for each gene during the analysis, to increase computational efficiency. Further details can be found in ‘Supplemental Methods—SNP-wise gene analysis’.

The MAGMA SNP-wise models can also be used to analyse raw genotype data, in which case the raw genotype data takes the place of the reference data and the SNP p-values are computed internally. Gene-set analysis based on these SNP-wise models proceeds in the same way as the gene-set analysis based on the multiple regression gene analysis model. The gene p-values resulting from the analysis are converted to Z-values in the same way to serve as input for the gene-set analysis. Similarly, the gene-gene correlation matrix R is obtained using the same formula as with the multiple regression model, but using the reference data to compute it.

Other features and implementation

A number of additional features has been implemented in MAGMA, more fully described in ‘Supplemental Methods—Extensions’. Gene analysis can be expanded with a gene-environment interaction component, which can subsequently be carried over to the gene-set analysis. Options for aggregation of rare variants and for fixed-effects meta-analysis for both gene and gene-set analysis are also available. Efficient SNP to gene annotation and a batch mode for parallel processing are provided to simplify the overall analysis process. MAGMA is distributed as a standalone application using a command-line interface. The C++ source code is also made available, under an open source license. MAGMA can be downloaded from http://ctglab.nl/software/magma.

Data

To evaluate the performance of MAGMA, the WTCCC Crohn’s Disease (CD) GWAS data [15] in conjunction with both WTCCC control samples was used. The data was cleaned according to the protocol described by Anderson [22], resulting in a sample of 1,694 cases and 2,917 controls with data for 403,227 SNPs. The European samples from the 1,000 Genomes data [23] and the HapMap 3 data [24] were used as reference data sets for the summary statistics gene analysis.

SNPs were annotated to genes based on dbSNP version 135 SNP locations and NCBI 37.3 gene definitions. For the main analyses only SNPs located between a gene’s transcription start and stop sites were annotated to that gene, yielding 13,172 protein-coding genes containing at least one SNP in the CD data. An additional annotation using a 10 kilobase window around each gene was made, yielding 16,970 genes, to determine the effect of using a window on relative performance. These two gene annotations were used for all analyses, to ensure that differences in default annotation settings did not cloud the comparison between tools. The 1,320 Canonical Pathways from the MSigDB database [16] were used for the gene-set analysis. The relatively large number of gene sets and the fact that the MSigDB Canonical Pathways are drawn from a number of different gene-set databases ensures a wide variety of gene sets, which should prevent the results from being too dependent on the choice of gene-set database.

Analysis of CD data—gene analysis

The results of the gene analyses of the CD data are summarized in Table 2, which shows the number of significant genes at a number of different p-value thresholds. Since the Type 1 error rates have been shown to be properly controlled these results can serve as a good indicator of the relative power of the different methods, and compared to simulation-based power estimates this has the advantage that no assumptions about the genetic causal model. From Table 2 it is clear that whereas the power of all the other methods is very similar, the MAGMA-main model shows a clear advantage over the rest. After Bonferroni correction, MAGMA-main found a total of 10 genome-wide significant genes, including the well-known CD genes NOD2, ATG16L1 and IL23R [25,26]. This also indicates that although MAGMA can perform analysis of summary statistics, raw data analysis should always be preferred if possible.

Specific implementation issues can be ruled out as the cause of the power difference since the PLINK and VEGAS analyses yield results highly similar to their matched MAGMA models (S9 Fig), and using the pruning option in PLINK also has little effect on the overall results. This means that the difference must be due to the difference in the methods and test-statistics themselves. Comparing the MAGMA implementations of these models in Fig 1, the mean χ ² and top χ ² approaches are shown to produce very similar p-values. Moreover, the plots reveal that the superior power of the MAGMA-main model does not arise from consistently lower gene p-values, but rather from a small set of genes with low p-values for MAGMA-main that are simply not picked up by the other approaches. This is likely to be related to the way LD between SNPs is handled, as that is one of the key differences between the multiple regression model of MAGMA-main and all the others. A post-hoc power simulation indeed indicates that multi-marker effects with weak marginals are the most probable explanation (see ‘Supplemental Methods—Simulation Studies’).

Open in a separate window

Fig 1

Comparison of gene analysis results for different test-statistics.

Gene -log10 p-values from the CD data gene analysis in MAGMA for three different gene test-statistics, comparing analyses using (A) the mean χ ² statistic with the top χ ² statistic, (B) the mean χ ² statistic and the PC regression model and (C) the top χ ² statistic and the PC regression model. P-values below 10^–8 are truncated to 10^–8 (grey points) to preserve the visibility of the other points.

To increase the generalizability of these findings, two variations on the CD analyses were performed for MAGMA-main, MAGMA-mean and MAGMA-top. First, the analyses were repeated with 10 principal components computed from the whole data set as covariates to correct for possible stratification. The results are shown in Table 2 and S10 Fig. There is shown to be only very limited stratification, and although the power does decrease somewhat MAGMA-main’s power advantage is maintained. The analyses were also repeated with the gene annotation extended to include a 10 kilobase window around each gene, with the comparison in S11 Fig showing a considerable impact on the results. However, although this suggests that the choice of window can strongly affect the results of a gene analysis Table 2 shows that the relative power stays the same, with MAGMA-main again maintaining its superior power.

Analysis of CD data—gene-set analysis

As with the gene analysis, the results of the CD analysis (Table 3 and Fig 2) can again serve as a gauge of the relative power of the different gene-set analysis methods. For the self-contained gene-set analysis this comparison is straightforward with MAGMA showing considerably more power than the two PLINK analyses. For the most part MAGMA’s power advantage can be explained by the difference in the underlying gene model, given the superior power of the PC regression model over the SNP-wise model used by PLINK shown before. Differences in how the genes are combined may also play a role however since, in contrast to PLINK, MAGMA weighs genes equally rather than by the number of SNPs in them and explicitly takes correlations between genes into account. Of note is also that PLINK-prune does considerably better than PLINK-avg, and that its p-values are somewhat more strongly correlated with those of the MAGMA analysis (Fig 2). An additional summary statistics analysis (MAGMA-pval-1K) on SNP p-values and using 1,000 Genomes reference data was also performed. This showed less power than PLINK even though it uses the same model at the gene level, suggesting that the difference is due to how the genes are aggregated to gene-sets. One of the key differences in this regard is that PLINK gives larger genes greater weight whereas MAGMA weighs them equally. As such a likely explanation is that the PLINK results are partially driven by a smaller number of large genes, though constructing the intermediate models to verify this is beyond the scope of this paper.

Open in a separate window

Fig 2

Comparison of self-contained gene-set analysis results.

Gene set—log10 p-values from the CD data self-contained gene-set analysis for MAGMA and PLINK. Panel (A) shows the PLINK-avg (no pruning) results compared with the MAGMA-main analysis, panel (B) the PLINK-prune results compared with the MAGMA-main analysis and (C) the two PLINK analyses compared to each other. P-values below 10^–8 are truncated to 10^–8 (grey points) to preserve the visibility of the other points.

The comparison of competitive methods is somewhat more complicated, due to the fact that ALIGATOR, INRICH and MAGENTA all use discretization using a p-value cut-off. This cut-off needs to be specified by the user and has no obvious default value, although for MAGENTA the 5^th percentile cut-off is suggested as the most optimal [12]. For ALIGATOR and INRICH the analysis was therefore performed at four different cut-offs (0.0001, 0.001, 0.005, 0.01), and for MAGENTA at two (5^th and 1^st percentile).

Of the four tools, only MAGMA and INRICH yield significant results after multiple testing correction (Tables ​(Tables33 and ​and4).4). As with the self-contained gene-set analysis, power for the MAGMA analysis is better when using raw data rather than SNP p-values as input, though both yield one significant gene set. For INRICH the results are strongly dependent on the SNP p-value cut-off used, with three significant gene sets at the 0.0001 cut-off but none at the higher ones, further emphasizing the problem of choosing the correct cut-off. It should also be noted that the p-values have not been corrected for the fact that the gene-sets have been analysed under four different thresholds, and thus might not fall below the significance threshold if they were.

Looking at the results in more detail (Fig 3) also suggests that the differences in results are not merely due to a difference in power. The concordance between methods is poor, with only MAGENTA and ALIGATOR showing a reasonable correlation in results. Moreover, there is considerable discordance between different p-values cut-offs for the same methods as well (Fig 4). This suggests that the different methods, or methods at different p-value cut-offs, are sensitive to distinctly different kinds of gene set associations. In particular, MAGMA and the other three methods at higher p-value cut-offs would be expected to respond best to gene-sets containing a larger number of somewhat associated genes. Conversely, at lower p-value cut-offs the latter three should become more sensitive to gene-sets containing a small number of more strongly associated genes. This is exemplified by the INRICH analysis. At the 0.0001 cut-off only quite strongly associated genes are counted as relevant, but as there are only 42 such genes overall the three gene sets (containing either 26 or 29 genes) become significant despite each containing only three relevant genes.

Open in a separate window

Fig 3

Comparison of competitive gene-set analysis results.

Gene set -log10 p-values from the CD data competitive gene-set analysis for MAGMA, ALIGATOR, INRICH and MAGENTA. Results for ALIGATOR and INRICH are shown for each for the SNP p-value cutoff that yielded the highest observed power (0.01 and 0.0001 respectively), MAGENTA at the advised 5^th percentile cutoff. P-values for gene sets not evaluated by one of the methods are shown in grey. The shown correlations are for the -log10 p-values for gene-sets evaluated by both methods.

Open in a separate window

Fig 4

Comparison of competitive gene-set analysis results at different SNP cut-offs.

Comparison of gene set -log10 p-values from the CD data competitive gene-set analysis at different SNP p-value cut-offs for ALIGATOR (top row), INRICH (middle row) and MAGENTA (bottom row). The highest cut-off on the horizontal axis is compared to each of the lower cut-offs. P-values for gene sets not evaluated at the lower cut-off are shown in grey. The shown correlations are for the -log10 p-values for gene-sets evaluated at both cut-offs. Horizontal and vertical grey dotted lines demarcate the p = 0.05 nominal significance threshold.

Aside from differences between methods, Table 3 also shows a clear difference between self-contained and competitive gene-set analysis. This is not a difference in power, but rather a difference of null hypothesis. Competitive tests attempt to correct for the baseline level of association present in the data and accordingly have a much more general null hypothesis. The impact of this difference in hypothesis can be illustrated by comparing the MAGMA self-contained and competitive analyses, since they are performed in the same framework. Whereas the self-contained analysis detects 39 gene sets that show association with the phenotype, the competitive analysis detects only one of those 39. For the remaining 38 gene sets, there is no evidence in the data that the associations in those gene sets are any stronger than would be expected by chance given the polygenic nature of CD. The gene-set that remains is the Regulation of AMPK via LKB1 (REACTOME) set. For two additional gene sets, Cell Adhesion Molecules (KEGG) and ECM-receptor Interaction (KEGG), the competitive p-value also drops below the significance threshold (Table 4 and S12 Fig) if the correction for gene size and gene density is turned off. This suggests that these gene sets do in fact contain significantly elevated levels of association, but that this is partially caused by confounding effects of the size and density of the genes they contain. Given the strength of the confounding effect it is evident that gene-set analyses should always be corrected for these and other potential confounders, to avoid false positive results. Full results for the analyses can be found in Table S5 in S2 File.