Comparisons of seven algorithms for pathway analysis using the WTCCC Crohn's Disease dataset

Simple single-marker tests used in genome-wide association studies (GWAS) have contributed to the discovery of many loci responsible for the variation observed in complex traits or disorders [1–3]; nevertheless, they are also criticized for their stringent significance threshold [4] and disregard of prior knowledge, which might lead to Type II errors, that is, not detecting real effects. Recently, a number of complementary approaches have been developed to prioritize susceptibility genes of complex diseases and increase power, of which meta-analysis, epistasis analysis and GWAS pathway analysis (GWASPA) are typical ones and already widely applied [5].

Biological pathways, which are actually series of actions among molecules in a cell that lead to a certain product or a change in a cell [6], are usually identified by experimental approaches [7, 8] and then revised with bioinformatical mining tools [9]. Shared across different organisms, pathways play an important role in metabolism, gene regulation and signal transduction [10]. With genetic or environmental perturbations, some normal pathways might become dysfunctional and then contribute to complex diseases [11]. Pathway analysis for complex disease, also known as gene-set analysis, originated with genomic expression studies [12]; one of the two classical methods in that field is over/under representation analysis by hypergeometric test, while the other is gene-set enrichment analysis (GSEA) using Kolmogorov-Smirnov-like test statistics [13, 14]. Since Wang et al. [15] first applied GSEA to GWAS data, more and more algorithms for performing genome-wide pathway analysis for SNP-chip datasets have emerged [16–19]. The PLINK set-based test utilizes average test statistics of groups of independent and/or truncated SNPs to provide a pathway-level test [17]; gene set ridge regression in association studies (GRASS) assesses joint association of pre-selected Eigen-SNPs for each gene in a candidate pathway with disease [18]; improved GSEA and association list Gene-ontology (GO) annotator (ALIGATOR) are two algorithms which utilize SNP-level test statistics or p-values in order to reduce computation cost over raw data-based methods [19, 20]. These different algorithms could be categorized according to type of input data (raw data or summary statistics), counts of analysis stages (one or two stage), or basic hypothetical tests (competitive or self-contained) [16]. The difference between GWAS pathway analysis and classical pathway analysis in genomic expression studies lies in whether a gene-based score is directly provided or not. Previous research on multi-allelic association tests or gene-based tests, which always produce gene-level scores indirectly, provide a good bridge for pathway analysis on individual SNPs [21–23]. Among them, GATES is a rapid and powerful procedure for getting gene-based statistics which sometimes serve as a prerequisite for further advanced analysis [22].

Though several methods have already been proposed, there is no consensus as to the best method for conducting a GWASPA, especially when the underlying causal mechanism for disease at the functional level is not yet clear [16, 24]. Results discrepancies among different methods on the same dataset might arise due to different mapping strategies from SNPs to pathways, or an algorithm's power of detecting susceptibility pathways, as well as incompatible pathway databases used [24]. Wang et al. [16] discussed the basic issues, main procedures and challenges involved in method development for GWASPA, but did not give a guideline of how to apply these methods from a practical perspective. Chen et al. [18] did compare the performance of different methods when evaluating the GRASS algorithm, but the simulation scenarios were designed for only one candidate pathway, and therefore need to be extended by examining both comprehensive simulated and real datasets on a genome-wide scale. Ballard et al. [25] illustrated the advantage of the random set method over the hypergeometric test for pathway analysis by analyzing three GWAS samples for Crohn's Disease (CD); however, they did not include any summary statistics-based algorithms which are becoming more and more popular [26].

As one of those complex diseases investigated at the earlier stage of GWAS [27], evidence for disease-causal variants, genes or even pathways related to CD are increasingly provided [3, 25, 27–29]. Meanwhile, the CD dataset from the Wellcome Trust Case Control Consortium (WTCCC) [30], which is openly available, has been repeatedly utilized and extensively explored by a variety of approaches [29, 31, 32]. This relatively abundant knowledge makes the WTCCC CD dataset a good testing sample for evaluating the performance of newly developed methods. Therefore, in this study, we implemented two summary statistics-based algorithms in an open-source tool named Knowledge-Based Mining System for Genome-Wide Genetic Studies (KGG) [22]; and then compared the performance of these two methods with another five existing methods using the WTCCC CD dataset, to evaluate the characteristics of the various methods, including Type 1 error, power and running time.

It has been suggested that multiple genes in immune system functional pathways, especially those containing different Interleukin factors, might be involved in causation of Crohn's Disease [28, 29]. Our findings of CD pathways across different algorithms also fall mostly into the category of immune response related pathways, though the number of significant pathways varied greatly for each individual algorithm. Likely, these differences are due to differences in four factors. The first is null hypothesis differences between self-contained tests and competitive tests [44, 45]. PLINK-Ave, PLINK-Max, GRASS and GATES-Simes, which are self-contained tests (see table 4), assume that a pathway does not contain any significant SNPs or genes. But GATES-Hyper, Aligator and GSEAforGWAS are all based on competitive tests (see table 4), which aims to test whether genes in one pathway are enriched with a greater number of associated SNPs. The second contributor to the differences observed is differences in test statistic construction. PLINK-AVE and GRASS both applies an "average" concept which combines evidence of selected SNPs or genes to identify pathways with more overall association signal, but PLINK-Max, GATES-Simes, GATES-Hyper, Aligator and GSEAforGWAS all use a "maximum" concept when computing pathway score from SNP scores directly or indirectly, relying on the most associated SNP in a pathway. It is likely that the underlying causal mechanism for CD involves more pathways covering moderate effect SNPs than simply a few pathways with top SNP hits only [21], implying that many more SNPs than we can identify individually contribute to disease. A third contributor to the differences observed is whether genes are treated as an intermediate bridge from SNPs and pathways in a two-stage approach. Though it has been suggested that two stage pathway analysis was more robust and powerful than one stage pathway analysis [46], our analyses don't support this, with the PLINK set-based test found to be superior to other gene-based algorithms. Finally, different strategies of handling LD structure between SNPs can also explain part of the variation in performance. Pruning SNPs in high LD (PLINK), choosing Eigen SNPs by principal component analysis (GRASS and KGG), or focusing only on the top SNP (Aligator and GSEAforGWAS), are all effective ways of reducing inflation caused by SNP dependency, but may perform differently.

Due to complexity caused by SNP, gene, and pathway relationships, re-sampling strategies like permutation are always utilized to derive an empirical distribution for test statistics to evaluate the significance of candidate pathways. However, it comes with a cost of time and memory. Usually two approaches are used to lighten the computational burden, either a revised adaptive or optimal permutation of raw genotype/phenotype data [46, 47] or a permutation of SNP or gene labels in summary statistics data [48]. We showed that GATES-Hyper and GATES-Simes implemented in KGG were much faster than those permutation-based algorithms since both of them used approximation for p-value distribution. Clearly, the advantages of GWASPA over genomic expression studies or traditional laboratory works (such as knock out mouse and cultured cells) [7, 8] are apparent in terms of time and economy and can be used to guide these further studies.

There are also some limitations of our study. While the WTCCC CD dataset is an ideal model GWAS, the performance of GWASPA algorithms may not be as good with other datasets, where studies are smaller and less is known regarding function at the outset. We recommend applying two or more algorithms from different categories (self-contained versus competitive, average versus maximum and one-stage versus two-stage) in practice, especially when we cannot readily distinguish the scenario in which various susceptibility variants confer moderate risk to disease versus the scenario in which a major effect variant in a pathway plays a dominant role in complex diseases [15, 20]. In addition, we set a fixed threshold for SNP or gene p-value truncation when including them for pathway analysis by PLINK-Ave, Aligator, and GATES-Hyper, but this cut-off is somewhat arbitrary and the choice made can influence results dramatically. Moreover, the FDR method was applied to correct for multiple testing and determine significance on pathway level, since almost no pathway p-values would survive Bonferroni adjustment [49]. Nevertheless, we think it was still reasonable for comparing performances of different algorithms, since we applied the same standards to all methods tested [24]. Lastly, none of these algorithms took differences in pathway structure stored in the original BioPAX or SBML format into consideration [50], but treated them as plain text containing independent gene symbols, when in fact any gene can be a member of multiple pathways. This drawback is handling by more advanced statistical methods like mixed-effect models and Bayesian networks, which facilitates modelling of gene-gene overlapping, interaction and correlation as well as net gene effect within the same pathway [51, 52].

GWAS pathway analysis, which prioritizes candidate pathways associated with complex disorders, could serve as an important complement to individual SNP analysis and gene-based analysis. Though all algorithms selected in this study do not have inflated Type I error rates, they vary greatly in terms of power and running time. The PLINK truncated set-based test was the most powerful, but the two summary statistics-based algorithms implemented in KGG were the fastest. However, raw data-based algorithms should be preferred for GWAS pathway analysis as long as computation capacity is available, since they preserve the intact data structure and tend to be more powerful than summary statistics-based algorithms. When underlying disease causal mechanism is ambiguous, which is common for complex diseases, it is worthwhile to apply two or more pathway analysis algorithms on the same GWAS dataset.

KGG is implemented using Java; therefore a Java Runtime Environment (JRE) is required to run KGG. Currently, installation of JRE and KGG are supported for Windows, Mac OS × and Linux. Three command files (run.win.bat, run.mac.sh and run.linux.sh) are provided for users to run KGG easily. A graphical user interface will automatically appear once initiating the command file. Documentation, source code, and precompiled binaries can be downloaded from http://bioinfo.hku.hk:13080/kggweb/home.htm.