Abstract

However, although various methods for gene and gene-set analysis currently exist, they generally suffer from a number of issues. Statistical power for most methods is strongly affected by linkage disequilibrium between markers, multi-marker associations are often hard to detect, and the reliance on permutation to compute p-values tends to make the analysis computationally very expensive. To address these issues we have developed MAGMA, a novel tool for gene and gene-set analysis. The gene analysis is based on a multiple regression model, to provide better statistical performance. The gene-set analysis is built as a separate layer around the gene analysis for additional flexibility. This gene-set analysis also uses a regression structure to allow generalization to analysis of continuous properties of genes and simultaneous analysis of multiple gene sets and other gene properties. Simulations and an analysis of Crohn’s Disease data are used to evaluate the performance of MAGMA and to compare it to a number of other gene and gene-set analysis tools. The results show that MAGMA has significantly more power than other tools for both the gene and the gene-set analysis, identifying more genes and gene sets associated with Crohn’s Disease while maintaining a correct type 1 error rate. Moreover, the MAGMA analysis of the Crohn’s Disease data was found to be considerably faster as well.

Author Summary

These methods can be used When the effects of individual markers is too weak to detect, Which is a common problem When studying polygenic traits. Moreover, gene-set analysis can provide additional insight into functional and biological mechanisms underlying the genetic component of a trait. Although a number of methods for gene and gene-set analysis are available however, they generally suffer from various statistical issues and can be very time-consuming to run. We have therefore developed a new method called MAGMA to address these issues, and have compared it to a number of existing tools. Our results show that MAGMA detects more associated genes and gene-sets than other methods, and is also considerably faster. The way the method is set up also makes it highly flexible. This makes it suitable as a basis for more general statistical analyses aimed at investigating more complex research questions.

Introduction

However, despite growing sample sizes, the genetic variants discovered by GWAS generally account for only a fraction of the total heritability of a phenotype [2,3]. More than anything, GWAS has shown that many phenotypes, such as height [4], schizophrenia [5] and BMI [6] are highly polygenic and influenced by thousands of genetic variants with small individual effects, requiring very large sample sizes to detect them.

In gene analysis, genetic marker data is aggregated to the level of whole genes, testing the joint association of all markers in the gene with the phenotype. Similarly, in gene-set analysis individual genes are aggregated to groups of genes sharing certain biological, functional or other characteristics. Such aggregation has the advantage of considerably reducing the number of tests that need to be performed, and makes it possible to detect effects consisting of multiple weaker associations that would otherwise be missed. Moreover, gene-set analysis can provide insight into the involvement of specific biological pathways or cellular functions in the genetic etiology of a phenotype. Gene-set analysis methods can be subdivided into self-contained and competitive analysis, with the self-con-tained type testing whether the gene set contains any association at all, and the competitive type testing whether the association in the gene set is greater than in other genes [7].

However, one concern with most existing methods is that they first summarize associations per marker before aggregating them to genes or gene sets. As demonstrated by Moskvina et al. this makes the statistical power strongly dependent on local linkage disequilibrium (LD) [14], and also reduces power to detect associations dependent on multiple markers.

These are often very computationally demanding, and since no parametric model is used it is often not made explicit which null hypothesis is being evaluated and what assumptions are made. This makes it more difficult to determine the properties of the analysis such as how the significance of a gene set relates to the significance of its constituent genes or whether the analysis corrects for a polygenic architecture. This complicates the interpretation of results and hampers comparison between results of different gene-set analysis methods.

MAGMA’s gene analysis uses a multiple regression approach to properly incorporate LD between markers and to detect multi-marker effects. The gene-set analysis is built as a distinct layer around this gene analysis, providing the flexibility to independently change and expand both the gene and the gene-set analysis. Both self-contained and competitive gene-set analyses are implemented using a gene-level regression model. This regression approach offers a generalized framework which can also analyse continuous gene properties such as gene expression levels as well as conditional analyses of gene sets and other gene properties, and which can be extended to allow joint and interaction analysis of multiple gene sets and other gene properties as well. More traditional gene analysis models are also implemented, for comparison and to provide analysis of SNP summary statistics.

Simulation studies were performed to verify type 1 error rates for MAGMA. The CD data set was then analysed using MAGMA and with five commonly used other tools for gene and gene-set analyses, specifically VEGAS [17] , PLINK [8], ALIGATOR [9], INRICH [10] and MAGENTA [12]. The results show that MAGMA has greater statistical power than the other methods, while also being considerably faster.

Materials and Methods

Model structure

In the first part a gene analysis is performed to quantify the degree of association each gene has with the phenotype. In addition the correlations between genes are estimated. These correlations reflect the LD between genes, and are needed in order to compensate for the dependencies between genes during the gene-set analysis. The gene p-values and gene correlation matrix are then used in the second part to perform the actual gene-set analysis. The advantage of decoupling these two parts of the analysis in this manner is that each can be changed independently from the other, simplifying the development of changes and extensions to either part of the model. Moreover, since the second part only uses the output from the first part and not the raw genotype data they do not need to be performed at the same time or place, making it much more straightforward to perform multiple gene-set analyses on the same data or to analyse multiple data sets across a large-scale collaboration.

Gene analysis

This model first projects the SNP matriX for a gene onto its principal components (PC), pruning away PCs with very small eigenval-ues, 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; the matriX of PCs, Y the phenotype and W an optional matriX of covariates the model can thus be written as Y = 050g —|— X; ag —|— Wflg —|— 5g, where the parameter vector 05g represents the genetic effect, fig the effect of the optional covariates, 060g the intercept and 6g the vector of residuals. The F-test uses the null-hypothesis H0: ag = 6 of no effect of gene g on the phenotype Y, conditional on all covariates.

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.

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 Meth-ods—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 pg computed with the gene analysis is converted to a Z—value 2g 2 (ID—1(1 — pg), where (ID—‘1L is the probit function. This yields a roughly normally distributed variable Z with elements zg 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 5 of the form Z5 2 flOT —|— 55, where Z5 is the subvector of Z corresponding to the genes in 5. Evaluating fig 2 0 against the alternative fig > 0 yields a self-contained test, since under the self-contained null hypothesis that none of the genes is associated with the phenotype zg 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 SS with elements 5g is then defined, with 5g 2 1 for each gene g in gene set 5 and 0 otherwise. Adding SS as a predictor of Z yields the model Z = flOST —l— Ssfls —l— 5. The parameter A 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 fis = 0 against the one-sided alternative A > 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 2 Q; —|— C1 Bl —|— Czflz —|— . . . —|— 5. The variables C1, C 2, . . ., 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 A 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. 5~MVN(6>, 021 However, due to LD, neighbouring genes will generally be correlated, violating this assumption. This issue can be addressed by using Generalized Least Squares ap-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 ‘Supple-mental Methods—Implementation Details’). Note that for the self-contained analysis, the sub-matriX R5 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 zg have a standard normal distribution, the error variance 02 can be set to 1.

Analysis of summary SNP statistics

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].

For the mean 12 statistic, a gene p-value is then obtained by using a known approximation of the sampling distribution [20,21]. For the top 12 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 12 statistic is computed for every gene. The empirical p-value for a gene is then computed as the proportion of permuted top 12 statistics for that gene that are higher than its observed top 12 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’.

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-val-ues 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

Gene analysis can be expanded with a gene-environ-ment 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

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.

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

The MAGMA gene analysis was performed on the raw CD data using the PC regression model

Gene analyses with VEGAS and PLINK were performed using the mean SNP statistic for VEGAS and both the mean SNP statistic (PLINK-avg) and the top SNP statistic (PLINK-top) for PLINK. Pruning in PLINK was turned off for these analyses. An additional PLINK analysis using the mean SNP statistic with pruning set to its default (PLINK-prune) was performed as well. To facilitate the comparison, several additional SNP-wise gene-set analyses were performed in MAGMA with test-statistics matching those of PLINK-avg, PLINK-top and VEGAS: mean [2 (MAGMA-mean) and top [2 (MAGMA-top) on the raw CD data to match the two PLINK analyses, and mean 22 using CD SNP p-values and with either HapMap reference data (MAGMA-pval) to match VEGAS or with 1,000 Genomes reference data (MAGMA-pval- 1K). The SNP summary statistics used for VEGAS and MAGMA-pval were computed using PLINK ‘—assoc’.

Several other analyses were performed for comparison: PLINK self-con-tained gene-set analysis without pruning (PLINK-avg) and with pruning (PLINK-prune), as well as ALIGATOR, INRICH and MAGENTA competitive gene-set analysis. PLINK operates on raw genotype data, Whereas all three competitive methods require only SNP p-values as input. No correction for stratification was used in any of the analyses except when explicitly specified. An overview of all analyses is given in Table 1.

Results

Type 1 error rates Simulation was used to assess the type 1 error rates, using permutations of the CD phenotype to obtain a global null distribution of no associated SNPs (see ‘Supplemental Methods—Simu-lation Studies’ for details). For the gene analysis, type 1 error rates were found to be controlled at the nominal level of 0.050 for the PC regression model, the summary statistics analysis model, as well as the SNP-wise models (Table S1 in 82 File).

For the competitive test an additional simulation using a polygenic null model was performed, with effects explaining a combined 50% of the phenotypic variance assigned to randomly selected SNPs. This polygenic type 1 error rate was also well controlled. The type 1 error rates for the self-contained analysis under the polygenic null model are also shown. These are considerably inflated because self-contained gene-set analysis by its definition is not designed to correct for polygenicity, illustrating the risk of performing self-contained analysis on polygenic phenotypes.

Analysis of CD data—gene analysis

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 signif1cant 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 952 and top [2 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 ‘Supplemen-tal Methods—Simulation Studies’).

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 810 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 811 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

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.

This cutoff needs to be specified by the user and has no obvious default value, although for MAGENTA the 5th percentile cutoff is suggested as the most optimal [12]. For ALIGATOR and INRICH the analysis was therefore performed at four different cutoffs (0.0001, 0.001, 0.005, 0.01), and for MAGENTA at two (5th and 13t percentile).

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 cutoff used, with three significant gene sets at the 0.0001 cutoff but none at the higher ones, further emphasizing the problem of choosing the correct cutoff. 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.

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 cutoffs for the same methods as well (Fig 4). This suggests that the different methods, or methods at different p-value cutoffs, are sensitive to distinctly different kinds of gene set associations. In particular, MAGMA and the other three methods at higher p-value cutoffs would be expected to respond best to gene-sets containing a larger number of somewhat associated genes. Conversely, at lower p-value cutoffs 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 cutoff 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.

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-con-tained 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 LKBI (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 812 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.

Computational performance

In terms of computational performance MAGMA is shown to have a considerable advantage over the other methods (Table 5) for both gene and gene-set analysis. The most marked difference is between MAGMA and PLINK, the only one of the alternative methods using raw data input. However, the raw data analysis in MAGMA outperforms the summary statistics methods as well. Although INRICH and ALIGATOR show comparable computation times at their lowest SNP p-value cutoff, the need to repeat the analysis at multiple cutoffs means the total analysis for both takes considerably longer.

Since the statistical tests used have known asymptotic sampling distributions the need for com-putationally demanding permutation or simulation schemes is avoided. Note however that the permutation-based SNP-wise analyses in MAGMA also show very reasonable computation times. These results demonstrate that, given efficient implementation, there is no computational reason to prefer analysis of summary statistics over raw data analysis, even when using permutation.

Discussion

Comparison with a number of other, frequently used methods shows that MAGMA has better power for gene analysis as well as for both self-contained and competitive gene-set analysis. An important factor in this is the multiple regression model used in the gene analysis, which is better able to incorporate the LD between SNPs than other methods. Because of its two-layer structure, this improvement in power at the gene-level subsequently carries over to the gene-set analysis.

This is primarily due to the choice of statistical model, which did not require the kind of computationally expensive permutation or sampling procedures used in the other methods. However, even the permutation-based SNP-wise models implemented in MAGMA outperformed their equivalents in other software and yielded very reasonable computation times. Although MAGMA showed better power than other tools for both the self-contained and competitive gene-set analysis, these comparisons also revealed considerable differences between the methods. This was most pronounced for the competitive gene-set analysis, with even results for individual methods showing significant variability based on the choice of cutoff. At present no comprehensive evaluation of the differences between existing gene-set analysis methods exists, leaving the causes and implications of these difference unclear. It is beyond the scope of this paper to perform such an evaluation, but the degree of discordance between most methods strongly suggests a need for future research in this direction. An additional caveat is that it is unknown to what extent the observed differences in power between methods may depend on the specific genetic architecture of Crohn’s diseases, and as such generalizing the results to other genetic architectures must be done with caution.

Because of the two-tiered structure of the gene-set analysis, alternative gene analysis models are straightforward to implement and are automatically available for use in the gene-set analysis. Similarly, the linear regression structure used to implement the gene-set analysis offers a high degree of extensibility. At present it enables analysis of continuous gene-level covariates as well as conditional analysis of gene-sets correcting for possible confounders, and the analysis of the CD data demonstrates that correction for confounders such as gene size and gene density is indeed strongly advised. The model is easily generalized to much more general gene-level linear regression models to allow for simultaneous analysis of multiple covariates and gene-sets, opening up a wide range of new testable hypotheses.

Supporting Information

Results for all Crohn’s Disease gene and gene-set analyses. (XLSX)

Empirical p-values were obtained for the CD data PC regression gene analysis by permutation of the F-statistic (A), in order to verify the accuracy of the asymptotic F-test p-values. An initial 100,000 permutations were computed for each gene. For genes with a very low initial empirical p-value (shown in blue and red) the number of permutations was increased to about 500 million to refine the empirical p-value. The dashed horizontal line indicates the lowest possible nonzero permutation p-value, genes with an empirical p-value of 0 are shown at half that minimum p-value in the plot (in red). The process was repeated using a subsample of the CD data skewed 4:1 towards cases (B) or controls (C); and with evenly divided subsamples of N = 1000 (D), N = 500 and N = 250. Only the initial 100,000 permutations were performed for these analyses, genes with an empirical p-value of 0 are again shown at half the minimum nonzero p-value (in blue). (TIFF)

Gene p-Values were computed using a logistic regression model to compare against the linear regression model used in MAGMA. P-Values were computed using either a Score test (A) or a Likelihood Ratio test (B). Because the Likelihood Ratio test appeared to have significantly more power than both the Score test and the MAGMA F-test, empirical p-Values for the Likelihood Ratio test were computed by generating up to 10,000 permutations of the Likelihood Ratio statistic. This was compared to the asymptotic Likelihood Ratio test p-Values (C), revealing a downward bias in the asymptotic p-Values. The empirical p-Values were then compared to the MAGMA F-test p-Values (D), which shows that the apparent power advantage of the Likelihood Ratio test in (B) was due to the bias in the

The pruning implemented in MAGMA was applied to the genes in the CD data at different levels of the prune factor f (default is 0.999), which reflects the proportion of the total variance in the raw genotype data that is retained after pruning. The original number of genotyped SNPs in each gene is plotted against the number of PCs retained after pruning. The regression slope gives an estimate of the average proportion of PCs to SNPs.

The PLINK—indep option was used to obtain an estimate of the number of independent SNPs at different R2 values. The number of PCs retained by MAGMA at different values of the pruning factor f is plotted against the number of independent SNPs at the R2 value that provided the closest match. (TIFF)

MAGMA needs to impute missing genotype values in order to run the multiple regression model, which is done by single imputation using flanking SNPs. To validate this procedure a subset of genes was selected from the CD data, and genotype values in those genes were set to be missing for a specified fraction of all the genotype values (up to 10%), and gene p-Values were then computed after using the imputation to fill in those missing values. Gene p-Values were also computed for the original full data. For each fraction, missing data was simulated 100 times for each gene, and the 5th (black) and 95th (blue) quantiles of the p-Values of each gene were computed and plotted against that gene’s full data p-Value.

Gene analysis was performed on the CD data, and a joint empirical distribution gene SSM values was generated using 4,611 permutations of the phenotype (since the sample size of the CD data is 4,611). The correlation matrix was then computed from this distribution. In addition, a correlation matrix for 13,172 uncorrelated genes was simulated by generating 4,611 permations for 13,172 genes and computing the correlation matrix. This provides the distribution of correlation coefficients that would be expected if the genes were uncorrelated. A QQ-plot of these expected correlation coefficients are plotted against the observed correlation coefficients in (A), showing a clear surplus of high positive correlations for the CD data genes. A QQ-plot using only correlations between genes more than 5 megabases apart (B) reveals that this is due to short-range correla-S7 Fig. Visualisation of the gene Z-statistic correlation matrix for chromosomes 5 and 6. Gene analysis was performed on the CD data, and a joint empirical distribution of the gene SSM values was generated using 4,611 permutations of the phenotype (since the sample size of the CD data is 4,611). The correlation matrix for chromosomes 5 and 6 was plotted, with individual pixels corresponding to a pair of genes and the color (from white to black) proportional to the absolute value of the correlation between those genes. Correlations with absolute value smaller than 0.05 are set to 0 to reduce noise. The yellow area corresponds to genes within 5 megabases of each other, corresponding to gene pairs for which MAGMA computes the correlations (correlations between more distant genes are assumed to be 0); the dashed lines indicate the boundary between the two chromosomes. (TIFF)

Summary statistics gene analysis of CD data SNP p-Values was performed using different reference data-sets, using the SNP-wise mean 22 model. This was compared to the same SNP-wise analysis performed on the raw CD genotype data. Grey points correspond to genes not covered by the reference dataset. The reference data-sets used are (A) the CD data itself, (B) the 1,000 Genomes European panel (97 missing genes), (C) the HapMap 3 European panel (375 missing genes) and (D) the HapMap 3 African panel (623 missing genes).

Gene —log10 p-Values from the CD data gene analysis for equivalent gene test-statistics implemented in different tools. The gene test-statistics used are (A) the mean 12 statistic in MAGMA and PLINK, (B) the top 12 statistic in MAGMA and PLINK, (C) the mean 12 statistic in MAGMA and VEGAS with analysis based on SNP p-Values and HapMap 3 reference data and (D) the mean 12 statistic in MAGMA on raw data and with analysis based on SNP p-Values and Hap-Map 3 reference data. (TIFF)

Gene —log10 p-Values from the CD data gene analysis for the three MAGMA gene analysis models with 10 PCs as covariates to correct for stratification, and without. P-Values below 10'8 are truncated to 10—8 (grey points) to preserve the visibility of the other points.

Gene —log10 p-Values from the CD data gene analysis for the three MAGMA gene analysis models with additional 10 kilobase window around the transcription start and stop sites, and without. Genes only present in the 10 kilobase window analyses are omitted. P-Values below 10'8 are truncated to 10'8 (grey points) to preserve the visibility of the other points.

Gene —log10 p-values from the CD data analyses. When the correction is turned on (the default setting), the gene-set effect is conditioned on gene size and gene density. Grey dashed lines represent the Bonferroni-corrected significance threshold. The effective size of the gene (number of PCs in the gene after pruning) is used as a measure of gene size, the ratio of effective size and total number of SNPs as a measure of gene density. The correction is achieved by entering gene size and gene density, as well as the log of both, as predictors in the generalized gene-set analysis model alongside the gene-set indicator variable. (TIFF)

Author Contributions

Performed the experiments: CAdL. Analyzed the data: CAdL. Wrote the paper: CAdL IMM TH DP.