GPA: A statistical approach to prioritizing GWAS results by integrating pleiotropy information and annotation data

Genome-wide association studies (GWAS) suggests that a complex disease is typically affected by many genetic variants with small or moderate effects. Identification of these risk variants remains to be a very challenging problem. Traditional approaches focusing on a single GWAS dataset alone ignore relevant information that could potentially improve our ability to detect these variants. We proposed a novel statistical approach, named GPA, to performing integrative analysis of multiple GWAS datasets and functional annotations. Hypothesis testing procedures were developed to facilitate statistical inference of pleiotropy and enrichment of functional annotation. We applied our approach to perform systematic analysis of five psychiatric disorders. Not only did GPA identify many weak signals missed by the original single phenotype analysis, but also revealed interesting genetic architectures of these disorders. We also applied GPA to the bladder cancer GWAS data with the ENCODE DNase-seq data from 125 cell lines and showed that GPA can detect cell lines that are more biologically relevant to the phenotype of interest.


Introduction
Hundreds of genome-wide association studies (GWAS) have been conducted to study the genetic bases of complex human traits/diseases. As of Dec., 2013, more than 11,000 single-nucleotide polymorphisms (SNPs) have been reported to be significantly associated with at least one complex traits/diseases (see the web resource of GWAS catalog [11] http://www.genome.gov/ gwastudies/). Despite of these successes of GWAS, these significantly associated SNPs can only explain a small portion of genetic contributions to complex traits/diseases [18]. For example, human height is a highly heritable trait whose heritability is estimated to be around 80%, i.e., 80% of variation in height within the same population can be attributed to genetic effects [39]. Based on large-scale GWAS, about 180 SNPs have been reported to be significantly associated with human height [1]. However, these loci together can only explain about 5-10% of variation in height [1,18,37]. This phenomenon is referred to as "the missing heritability" [16,17,18].
Finding the missing heritability has drawn much attention from worldwide researchers. The role of common variants (minor allele frequency (MAF) ≥ 0.01) has been shown to be critical in explaining the phenotypic variance. Instead of only using genome-wide significant SNPs, Yang et al. [41] reported that, by using all genotyped common SNPs, 45% of the variance for human height can be explained. This result suggests that a large proportion of the heritability is not actually missing: given the limited sample size, many individual effects of genetic markers are too weak to pass the genome-wide significance, and thus those variants remain undiscovered. So far, people have found similar genetic architectures for many other complex traits/diseases [38], such as metabolic syndrome traits [36] and psychiatric disorders [5,12], i.e., the phenotype is affected by many genetic variants with small or modest effects, which is usually referred to as "polygenicity". The polygenicity of complex traits is further supported by recent GWAS with larger sample sizes, in which more associated common SNPs with moderate effects have been identified (e.g., [20]). Clearly, the emerging polygenic genetic architecture imposes a great challenge of identifying risk genetic variants: a larger sample size is required to identify genetic variants with smaller effect sizes. However, sample recruitment may be expensive and timeconsuming. Alternatively, integrative analysis of genomic data could be a promising direction, including combining GWAS data of multiple genetically related phenotypes and incorporating relevant biological information.
The last couple of years have seen an emerging recognition of "pleiotropy", i.e. the sharing of genetic factors, between human complex traits/diseases. For example, a systematic analysis of the NHGRI catalog (http://www.genome.gov/gwastudies/) of published GWAS showed that 16.9% of the reported genes and 4.6% of the reported SNPs are associated with multiple traits/diseases [30]. Through a "pleiotropic enrichment" method, Andreassen et al showed that it is possible to improve the power to detect schizophrenia-associated genetic variants by utilizing the pleiotropy between schizophrenia (SCZ) and cardiovascular-disease [2]. A more recent study identified four significant loci (p-value < 5 × 10 −8 ) to be pleiotropic by analyzing GWAS data of 33,332 cases and 27,888 controls for five psychiatric disorders [6]. Further analysis suggested very significant genetic correlation between schizophrenia and bipolar disorder (0.68 ± 0.04 s.e.) [5]. The wide-existing pleiotropy has also been demonstrated among several other types of traits, for example, metabolic syndrome traits [36] and cancers [26].
An increasing number of studies also suggest that functionally annotated SNPs are generally more biologically important and henceforth more likely to be associated with complex traits/diseases. To name a few, Schork et al. [27] demonstrated a consistent pattern of enrichment of GWAS signals among functionally annotated SNPs, using GWAS data of different traits/diseases (e.g., Crohn's disease and SCZ). Yang et al. [43] showed that SNPs in the genic region could explain more variance of height and body mass index (BMI) than SNPs in the intergenic region. Nicolae et al. [22] suggested that SNPs associated with complex traits/disease were more likely to be expression quantitative trait loci (eQTL). In addition, public availability of the vast amount of functional annotation data also provides unprecedented opportunities to investigate the enrichment of GWAS signals among these various types of functional annotations. For example, recently, the Encyclopedia of DNA Elements (ENCODE) Consortium has generated vast amounts of experimental data on gene expression (RNA-seq), DNA methylation status (RRBS-seq), chromatin modifications (ChIP-seq), chromatin accessibility (DNase-seq and FAIRE-seq), transcription factor (TF) binding sites (ChIP-seq), and long-range chromatin interactions (ChIA-PET, Hi-C, and 5C). As of September 2012, more than 1,600 data sets from 147 cell lines have been produced to annotate human genome, including 2.89 million unique, non-overlapping DNase I hypersensitivity sites (DHSs) in 125 cell lines using DNase-seq and 630K binding regions of 119 DNA-binding proteins in 72 cell lines using ChIP-seq, among many [34]. The ENCODE Project Consortium [34] examined 4,492 risk-associated SNPs from NHGRI GWAS catalog and found that 12% of them overlap with TF binding regions and 34% overlap with DHSs.
The widely existing pleiotropy between complex traits/diseases and increasingly rich functional annotation data are calling for novel statistical methods to effectively perform joint analysis of multiple GWAS data sets and functional annotation data. Statistical methods to investigate pleiotropy have been actively researched (reviewed in [32] and [29]), for example, using linear mixed models [14,15]. However, these methods do not allow to utilize functional annotation data for prioritization of GWAS results. On the other hand, various statistical methods which make use of functionally annotated SNPs have been proposed in recent years (reviewed in [8], [4], and [40]). For example, GSEA [33] identifies potentially important pathways in which target genes of risk-associated SNPs are involved while RegulomeDB [3] allows nucleotide-level annotations of risk-associated SNPs, especially for those located in non-coding regions. Stratified False Discovery Rate methods were applied to incorporate annotation into GWAS data analysis [27]. However, these methods are designed for the analysis of single phenotype and hence, these methods do not utilize functional annotation data fully efficiently for the genetic variants shared by multiple phenotypes. In short, well established statistical framework for the integration of functional annotation data to the joint analysis of genetically correlated GWAS data is still missing in the current literature.
In order to meet the emerging needs described above, we propose a unified statistical framework, named GPA 1 , to prioritize GWAS result based on the pleitropy and annotation information. GPA also provides statistically rigorous and easily interpretable inference tools for the analysis of pleiotropy and the enrichment of functional annotation. This article is organized as follows. In Section 2, we first describe our GPA probabilistic model and its statistical inference procedures. In Section 3, we investigate the properties of GPA using extensive simulation studies and illustrate the versatility and utility of GPA with the analysis of real data. Specifically, we apply GPA to the five psychiatric disorder GWAS data with central nervous system gene expression data and show that GPA can accurately identify pleiotropy structure among these diseases. We further apply GPA to the bladder cancer GWAS data with the ENCODE DNaseseq data from 125 cell lines and show that GPA can detect cell lines that are biologically more relevant to the phenotype of interest. In Section 4, we conclude this paper with the discussion of related issues.

GPA probabilistic model
Throughout this paper, we shall use j to index SNPs, k to index GWAS data sets, d to index the annotation data sets. Let's first consider the simplest case that we only have summary statistics (p-values) from only one GWAS data set, then extend our model to handle multiple GWAS data sets and annotation data. Suppose we have performed hypothesis testing of genome-wide SNPs and obtained their p-values: where M is the number of SNPs. Consider the "two-groups model" [10], i.e., the obtained p-values are assumed to come from the mixture of null and non-null, with probability π 0 and π 1 = 1 − π 0 , respectively. Let Z j = [Z j0 , Z j1 ] be the latent variables indicating whether the j-th SNP is null or non-null, where Z j0 ∈ {0, 1}, Z j1 ∈ {0, 1}, and Z j0 + Z j1 = 1 (A SNP can only be either null or non-null). Z j0 = 1 means un-associated (null) and Z j1 = 1 means associated (non-null). Then we have the following two-groups model: where the p-values from the null group are from the Uniform distribution on [0,1], denoted as U[0, 1], and the p-values from the non-null group are from the Beta distribution with parameters (α, 1), where 0 < α < 1. We put the constraint 0 < α < 1 to model that smaller p-value is more likely than larger p-value when it is from the non-null group [23]. To incorporate information from functional annotation data, we extend the basic model as follows. Suppose we have collected information from D functional annotation sources in the annotation matrix: A ∈ R M ×D , where A j,d ∈ {0, 1} indicates whether the j-th SNP is annotated in the d-th functional annotation source. For example, when there are two annotation sources: eQTL data and DNase I hypersensitivity sites (DHS) data, then A is an M × 2 matrix. If the j-th SNP is an eQTL, then A j1 = 1, otherwise A j1 = 0; if it is located in a DHS, then A j2 = 1, otherwise A j2 = 0. Now we model the relationship between Z j and A jd as Clearly, q d0 = Pr(A jd = 1|Z j0 = 1) can be interpreted as the proportion of null SNPs in the d-th annotation file, and q d1 = Pr(A jd = 1|Z j1 = 1) corresponds to the proportion of non-null SNPs in the d-th annotation file. Therefore, q d1 >> q d0 means that there exists enrichment in the d-th annotation file. The statistical inference about enrichment of annotation data will be discussed in details in Section 2.2.2. Now we extend the above model to handle multiple GWAS data sets. To keep the notation uncluttered, we present the model for the case of two GWAS data but the generalization to more than two GWAS data is straightforward. In fact, our GPA model and its algorithm are not limited by the number of GWAS data. Suppose we have p-values from two genetically related GWAS: p-value from GWAS1 :P 11 , P 21 , . . . , P j1 , . . . , P M 1 .
Let P ∈ R M ×2 be the matrix collecting all the p-values, where P jk denotes the p-value of the j-th SNP in the k-th GWAS. Similarly, we introduce latent variables Z j = [Z j00 , Z j10 , Z j01 , Z j11 ] indicating the association between the j-th SNP and the two phenotypes: Z j00 = 1 means the j-th SNP is associated with neither of them, Z j,10 = 1 means it is only associated with the first one, Z j01 = 1 means it is associated with the second one, and Z j11 = 1 means it is associated with both. The two-groups model (2) is extended to the following "four-groups model": where 0 < α k < 1, k = 1, 2. When the genetic basis of the two phenotypes are independent of each other (i.e., no pleiotropy), then we have π 11 = (π 10 + π 11 )(π 01 + π 11 ) by expectation. Therefore, the difference between π 11 and (π 10 + π 11 )(π 01 + π 11 ) can be used to characterize pleiotropy. Statistical inference on pleiotropy is given in Section 2.2.2.
To incorporate annotation information into multiple GWAS model (5), similarly, we model the relationship between Z j and A jd as where q d00 is the probability of a null SNP being annotated, q d10 is the probability of the first phenotype associated-SNP being annotated, q d01 is the probability of the second phenotype associated-SNP being annotated and q d11 is the probability of jointly associated-SNP being annotated. Assuming the independence of SNP markers, the joint distribution Pr(P, A) can be written as where P j and A j are the j-th row of P and A; the second equation holds by assuming the independence between P j and A j , conditional on Z jl ; the third equation holds by further assuming the independence between A jd and A jd for d = d , conditional on Z jl . Parameters of the GPA probabilistic model are estimated using the Expectation-Maximization (EM) algorithm [7], which turns out to be computationally efficient because we have explicit solutions for estimation of all the parameters in the M-step. Standard errors for parameter estimates are estimated using the empirical observed information matrix [19]. Note that in the GPA model, the sample size for estimating the empirical observed information matrix corresponds to the number of SNPs and as a result, we have a very large sample size (∼ 10 6 ) to accurately estimate standard errors. More details of the EM algorithm and the estimation of standard errors are provided in Appendix.

Statistical inference 2.2.1 False discovery rate
After we estimate parameters for the GPA model, SNPs can be prioritized based on their local false discovery rates [9]. Note that separate analysis of single GWAS data based on their summary statistics is equivalent to the analysis of single GWAS data using GPA without any annotation data. Hence, separate analysis of single GWAS data can be considered as a special case of our GPA framework.
To present the local false discovery rate based on our GPA model, we begin with the simplest case: single GWAS without annotation data, In this case, there are only two groups: null and non-null. The false discovery rate is defined as the probability that j-th SNP belongs to the null group given its p-value, i.e., fdr(P j ) = Pr(Z j0 = 1|P j ). (8) For joint analysis of two GWAS data sets, we are interested in the local false discovery rate of the j-th SNP, if it is claimed to be associated with the first phenotype and the second one, i.e., Similarly, when annotation data are available, the false discovery rates can be calculated as Then, we use the direct posterior probability approach [21] to control global false discovery rates.
More details for the estimation of false discovery rates are provided in Appendix.

Hypothesis testing of annotation enrichment and pleiotropy
Consider a given annotation file, we are particularly interested in the significance of its enrichment. First, we propose to the likelihood ratio test (LRT) to assess the significance. To keep the notation simple, we drop the index of annotation files. Specifically, the significance of enrichment of an annotation file for single-GWAS data analysis can be assessed by the following hypothesis: The LRT statistics is given as follows: is the parameter estimates obtained under H 0 (Here the superscript A indicates the Annotation enrichment test). Note thatΘ (A) 0 can be easily obtained by running the GPA algorithm without incorporating the annotation data. Under the null, the test statistic −2 log λ (A) asymptotically follows the χ 2 distribution with degree of freedom df = 1 [28]. We reject For joint analysis of two GWAS with annotation data, hypotheses (11) become as Under the null, the test statistics asymptotically follows χ 2 distribution with df = 3. Similarly, the test of annotation enrichment can be extended to handle K GWAS. In this case, the test statistics asymptotically follows χ 2 distribution with df = 2 K − 1 under the null. Now we consider testing pleiotropy between two GWAS. When there is no pleiotropy, i.e., the signals from the two GWAS are independent of each other, testing pleiotropy can be formulated by evaluating the following hypothesis: where π 1 * = π 10 + π 11 , π * 1 = π 01 + π 11 . The LRT statistics is constructed as follows: is the parameter estimates obtained under H 0 (Here the superscript P indicates the Pleiotropy test). The test statistics (−2 log λ (P ) ) asymptotically follows χ 2 distribution with df = 1 under the null.

Simulation study
We first evaluated performance of GPA for parameter estimation and hypothesis testing and compared them with some related methods using simulation studies. We followed the classical liability threshold model [12] to simulate case-control GWAS data for two genetically correlated diseases. For each disease, we first simulated a large cohort of individuals with genotypes of M independent SNPs. The MAFs of these SNPs were drawn uniformly from [0.05, 0.5]. Then we randomly chose m SNPs to be causal SNPs. The per-minor-allele effect of each causal SNP was drawn from a normal distribution with zero-mean and variance of where h 2 is the desired level of variance explained by all SNPs on the liability scale and f i is the MAF of the corresponding causal SNP. We also simulated the the environmental effect on the liability scale for each individual from a standard normal distribution (zero mean and unit variance). The total liability for each individual was then obtained by adding up all the genetic effects and the environmental effect. Given a desired disease prevalence B, individuals with liabilities greater than the 1 − B quantile were classified as cases and others were classified as controls. Then equal numbers of cases and controls were drawn from the cohort as a GWAS data set. When simulating two diseases simultaneously, we simulated two disjoint cohorts with the same set of SNPs. To reflect the pleiotropy effects between the two diseases, m causal SNPs (m ≤ m) were chosen to be shared by the two diseases. The annotation status of each causal and non-causal SNP was simulated from a Bernoulli distribution with probability of q 1 and q 0 respectively. In our simulation study, the sample size of each data set, N , was set at 1000, 2000, 5000 or 10000. The number of causal SNPs m was the same for the two diseases and was set at 500, 1000 or 2000. We varied the proportion of shared causal SNPs between the two diseases, γ, from m/M to 1. Note that γ = m/M corresponds to zero pleiotropy. The disease prevalence, B, was fixed at 0.1 and the variance explained by all the SNPs, h 2 , was fixed at 0.6 for the two diseases. Here q 1 and q 0 were fixed at 0.4 and 0.1, respectively.
We first evaluated SNP prioritization accuracy of GPA. Specifically, after the two data sets were simulated, we obtained the p-value for each SNP in each disease using a one degree of freedom χ 2 test. Then we analyzed the simulated data using our GPA method in the following four modes: 1. analyzing the two diseases separately without the annotation data; 2. analyzing the two diseases with the annotation data; 3. analyzing the two diseases jointly without the annotation data; 4. analyzing the two diseases jointly with the annotation data. The simulation was repeated 50 times for each scenario. In each mode, we compared the order of the local false discovery rates obtained using GPA against the actual causality status of the SNPs to Because all the simulation parameters are the same for the two diseases, only the results for the first disease is shown. The results suggest that both the annotation information and the pleiotropy between the two diseases improved the prioritization accuracy. In particular, as the fraction of shared causal SNPs increased, the prioritization accuracy also increased. Given the local false discovery rates obtained using GPA, we controlled the global false discovery rate at 0.2 and calculated the average power for all the true causal SNPs. The results for N = 5000 and m = 1000 are shown in Figure 1 (right panel) and are similar with results of the AUC (results for other scenarios are in the supplementary material, Figure S1-S11). We also checked the actual false discovery rates and found that on average the false discovery rate was indeed by controlled at 0.2 in spite of occasional slight conservativeness ( Figure 2, Figure S12-S22).
In the simulations, we also found that GPA provided satisfactory estimate of q, the probability of being annotated for a certain group of SNPs, as long as there are enough number of SNPs in that group (Figure 3, Figure S23-S33). However, we note that the proportion of causal SNPs always tends to be under-estimated ( Figure 4, Figure 5, Figure S34-S55). The underestimation may be due to the limited sample size and the model mismatch between GPA and the random-effects model. As a comparison, we also tried the "conditional FDR" approach proposed by Andreassen et al. [2] to prioritize the SNPs in our simulations. The comparison between this approach and GPA at N = 5000 and m = 1000 is shown in Figure 6. GPA significantly outperformed the conditional FDR approach in terms of the SNP prioritization accuracy. More importantly, when no pleiotropy exists, the conditional FDR approach showed worse accuracy than single-GWAS analysis using the standard FDR approach, whereas GPA achieved comparable accuracy with    The upper panel shows the estimated π 0 * , π 1 * , π * 0 , π * 1 with annotation, and the lower panel shows the estimated π 0 * , π 1 * and π * 0 , π * 1 without annotation. The red lines represent the true values.  Figure 5: The estimated π of joint analysis (two GWAS) at N = 5000, m = 1000. The upper panel shows the estimated π 00 , π 10 , π 01 , π 11 with annotation, and the lower panel shows the estimated π 00 , π 10 , π 01 , π 11 without annotation. The red lines represent the true values. single-GWAS analysis in this scenario. This suggests that GPA was able take advantage of pleiotropy while it sacrifices much less statistical power than the conditional FDR approach when pleiotropy does not exist. Next, we evaluated the power and type I error of GPA for hypothesis testing on the significance of annotation enrichment in causal SNPs. Another popular method that does a similar job is the Gene Set Enrichment Analysis (GSEA) [33]. Although GSEA typically is used for analysis of gene expression data, its input can be a list of p-values obtained from any type of sources. Therefore we implemented the GSEA method to test the enrichment the p-values of a set of SNPs being annotated and compared it with GPA. We followed the previous simulation scheme and simulated one GWAS data set with M = 20000, N varying from 2000 to 10000, m varying from 500 to 2000. Here q 0 was fixed at 0.1 and q 1 was varied from 0.1 to 0.5. We performed tests at the significance level of 0.05. Type I error rate was evaluated at q 1 = 0.1 and power was evaluated in other cases. The result for n = 1000 was shown in Figure 7. We observed that in general, GPA provided much higher power than GSEA while both methods appropriately controlled the type I error rate.
Lastly, we evaluated the power and type I error rate of GPA for the test of pleiotropy in our simulations. The simulation parameters were the same as in the previous simulations. Power was evaluated at γ = 0.25, 0.5, 0.75 and 1. The type I error rate was evaluated at γ = m/M . As shown in Figure 8, the power increases as m decreases and as N and γ increases, whereas the type I error rate was appropriately controlled in all cases.    Table 1: Single-GWAS analysis of five psychiatric disorders using the CNS gene set as the annotation data. Hereα is the estimate of α parameter of Beta distribution (2),π 0 andπ 1 are the estimated proportion of null-SNPs and non-null-SNPs defined in (2),q 0 andq 1 are the estimated proportion of null and non-null SNPs in the CNS gene set. Enrichment foldq 1 /q 0 and p-value given by hypothesis testing of enrichment in annotation data are provided in last two columns. The values in the brackets are standard errors of the estimates.

GWAS of five psychiatric disorders
We applied our GPA model to analyze the five psychiatric disorders [5,6]: attention deficithyperactivity disorder (ADHD), autism spectrum disorder (ASD), bipolar disorder (BPD), major depressive disorder (MDD) and schizophrenia (SCZ). Detailed information about these data sets is provided in [5,6]. The summary statistics of the five psychiatric disorders were downloaded from the section of cross-disorder at the website of Psychiatric Genomics Consortium (PGC) https://pgc.unc.edu/Sharing.php. The p-values were available for 1,230,535 SNPs in ADHD, 1,245,864 SNPs in ASD, 1,233,533 SNPs in BPD, 1,232,794 SNPs in MDD, and 1,237,959 SNPs in SCZ, respectively. We took the intersection of those SNPs and obtained their p-values. Finally, we had a p-value matrix P ∈ R 1,219,805×5 for the five psychiatric disorders. First, we performed single-GWAS data analysis using genes preferentially expressed in the central nervous system (CNS) [12,24] as the annotation data. Specifically, we generated the annotation vector A ∈ R 1,219,805 as follows: The entries in A corresponding to SNPs exactly in the genes from the CNS set or within 50-kb boundaries of those genes were set to be 1. Consequently, 21.9% of the genome-wide SNP markers were annotated. The analysis results of these five psychiatric disorders were given in Table 1 [5], suggesting about 1.6% and 1.5% fold enrichment in BPD and SCZ, respectively. The minor difference of fold enrichment estimation between GPA and LMM could be understandable: First, GPA only used summary statistics while LMM used both phenotype and genotype data. Second, the mathematical definition of fold enrichment is different between GPA and LMM: GPA used the ratio betweenq 1 andq 0 , while LMM used the ratio between the proportion of the variance explained by SNPs in the CNS set and the proportion of the CNS set in entire genome. Furthermore, we evaluated the significance of enrichment of the CNS set by hypothesis testing. As given in Table 1, enrichment of the CNS gene set was strong in BPD and SCZ, moderate in ASD and MDD, and nonsignificant in ADHD.
Next, we applied GPA to study pairwise pleiotropy of these five psychiatric disorders without  using annotation data. As shown in Table 2 , and non-significant for other pairs. Our analysis result agrees with the most of pleiotropy results reported in [5] and the disagreement mainly came from the joint analysis between ADHD and other disorders. The pleiotropy between ADHD and MDD was reported to be moderate, while GPA did not detect this moderate effect. From single GWAS analysis of ADHD, given in Table 1, the estimated parameters (π 1 =0.009 (s.e. 0.006),α=0.694) indicates that its GWAS signal was very weak. For MDD, the estimated parameterα=0.837 also indicates the weak marginal signal of MDD. Consequently, the marginal GWAS signals of ADHD and MDD were too weak to allow GPA to detect the pleiotropy effect between them.
Since the data analysis performed in [5] used genotype data, the bivariate linear mixed model could still have enough power to detect the moderate genetical correlation between ADHD and MDD.
We further applied GPA to study all pairs of disorders using the CNS gene set as the annotation data. The estimatedq l (l ∈ {00, 10, 01, 11}) are given in Table 3 andπ l remained almost the same as without annotation data. The p-values of hypothesis test (12) are also provided in the last column of Table 3. The p-value should be interpreted with caution: as shown in Table 1, the CNS gene set is enriched in all these disorders except ADHD. Hence, the significant p-values listed in Table 3 may be simply due to the combinatorial effects. On the other hand, the ratio betweenq 11 andq 00 could be more interesting. Take BPD-SCZ as an example. The ratio betweenq 11 andq 00 is 1.503 (s.e. 0.025), which suggests that enrichment of the CNS set for the BDP-SCZ shared risk variants was even stronger than that for BPD-only (1.467 (s.e. 0.033)) or SCZ-only (1.391 (s.e. 0.022)).
We also compared the results given by four different analysis approaches: single-GWAS analysis with or without annotation, two-GWAS joint analysis with or without annotation data. The manhattan plots are shown in Figure 9. For single-GWAS analysis without annotation, GPA identified 13 SNPs and 391 SNPs with f dr < 0.05 for BPD and SCZ, respectively. By using the CNS set as annotation, GPA was able to identify 14 and 409 SNPs for BPD and SCZ, respectively. For joint analysis without annotation, the number of identified SNPs dramatically increased to 383 and 821 for BPD and SCZ, respectively. By using the CNS set  Table 3: GPA results for all pairs of the five psychiatric disorders with the CNS set as annotation. The estimatedπ l almost remain the same as Table 2 andq l are shown in the table. The values in the brackets are standard errors of the estimates. The p-values of hypothesis testing (12) are provided in the last column.
as annotation, the number of identified SNPs further increased to 385 and 837 for BPD and SCZ, respectively. We investigated the result of BPD more in detail to evaluate the power of GPA in identification of functionally important SNPs. For single-GWAS analysis of BPD, GPA was able to identify SNPs locating the ANK3 gene. By using annotation data, the CACNA1C gene which encodes an alpha-1 subunit of a voltage-dependent calcium channel was identified by GPA. After incorporating pleiotropy information between SCZ and BPD, more functionally relevant genes, such as PBRM1, C6orf136, DPCR1, SYNE1, can be further identified by GPA. For instance, SYNE1 encodes the synaptic nuclear envelope protein 1, and provides instructions for making a protein called Syne-1 that is found in many tissues and especially critical in the brain. The Syne-1 protein is active (expressed) in Purkinje cells, which are located in the cerebellum and are involved in chemical signaling between nerve cells (neurons). Mutations in the SYNE1 gene have been found to cause autosomal recessive cerebellar ataxia type 1 (ARCA1) and SYNE1 has recently been implicated as a susceptibility gene for BPD in a large collaborative GWAS study [31]. Clearly, these results indicate that the statistical power to identify associated SNPs increased a lot by making use of pleiotropy and functional annotation (in this real data example, pleiotropy played a more important role than functional annotation). We also tried to apply GPA for joint analysis of BDP, SCZ and MDD to further explore their joint genetic architecture. In this scenario, GPA used eight-group to model all 2 3 states of a SNP, which is a nature extension of four-group model (5). However, we found that the estimated parameters obtained from this model were not reliable because of their large standard errors. Also, the interpretation of the estimated parameters became more complicated. Therefore, we mainly focused on two-GWAS analysis in this paper.
Regarding to computational time, the GPA algorithm takes less than 20 minutes to analyze typical GWAS data sets. The speed of convergence depends on the strength of the GWAS signals. For example, it took about 7 mins and 3 mins to analyze ADHD and SCZ, repectively, as SCZ has a stronger GWAS signal than ADHD. For joint analysis of BPD and SCZ, it took about 20 mins. All timings were carried out on a desktop with 3.0 GHz CPU with 16G memory.

Bladder cancer GWAS and ENCODE annotation data
In molecular biology, DNase I hypersensitive sites (DHSs) are regions where DNA degradation by enzymes like DNase I raise more frequently. As a result, DHSs can mark active transcription regions across genome and these patterns are known to be tissue or cell specific. The ENCODE project analyzed the DHSs in 125 human cell lines with the intention of cataloging human regulatory DNA [35]. In this section, we applied GPA to assess how bladder cancer [25] risk associated SNPs are enriched in DHSs region across these 125 human cell lines.
We downloaded genotype data for bladder cancer from dbGaP (NCI Cancer Genetic Markers of Susceptibility (CGEMS) project; accession number phs000346.v1.p1). We used the samples genotyped from both Illumina 1M chip and 610K chip for our analysis. For quality control, we removed SNPs with a missing rate > 0.01. We checked Hardy-Weinberg Equilibrium and excluded SNPs with p-value < 0.001. SNPs with minor allele frequency (MAF) < 5% were also removed. Finally, 490,614 SNPs from 3,631 cases and 3,356 controls of European descent were used in the analysis. Then, we analyzed this bladder cancer data with logistic regression by assuming additive model and obtained p-values. We also downloaded the uniform peak files for DHS in 125 cell lines from the ENCODE database (http: //genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeAwgDnaseUniform). Note that the DHS for these 125 cell lines were identified with a uniform analysis workflow by the EN-CODE Consortium and this facilitates fair and unbiased comparison among cell lines as annotation for our GPA model.
We applied GPA to analyze the bladder cancer GWAS data set with one annotation data at a time, and performed hypothesis testing to assess the significance of enrichment. The result is shown in the left panel of Figure 10. Under significance level α = 0.05 after Bonferroni correction, annotations from 19 cell lines were significantly enriched for bladder cancer risk associated SNPs. Most of these cell lines were derived from lymphocytes from normal blood (e.g., T cells CD4+ Th0 adult, Monocytes CD14+ RO01746), while some cell lines came from cancer patients (e.g., Gliobla and HeLa-S3). The above results demonstrates that the functional roles of bladder cancer risk variants might be involved in disturbance of immune system or carcinoma pathways. This also implies that GPA may be an effective way to explore functional role of GWAS hits by testing enrichment on phenotype-related annotations or user-specified annotations.
We also compared GPA with the LMM-based approach [13,43]. Specifically, we considered the following genome-partitioning linear mixed model: where X are covariates (the first five principal components from genotype data), and W 1 and W 2 are sets of SNPs overlapping DHS in each cell line and the remaining SNPs, respectively. We denote number of SNPs in W 1 and W 2 as M 1 and M 2 , respectively. Median number of SNPs that overlap DHS in each cell line is about 60K and 90% of cell lines have number of DHS ranging between 40K and 80K. In order to take into account such variation in number of DHS among cell lines, we defined a scaled version of the proportion of phenotype variance explained by SNPs overlapping DHS in each cell line as where is the proportion of the explained variance and M M 1 is the scaling factor. The right panel of Figure 10 shows that (−log 10 )-transformed p-value of the GPA annotation enrichment test is well linearly related to v 1 . This indicates that our GPA model captures enrichment of annotation almost as accurately as LMM even without the original genotype data, which implies wide applicability of our GPA model.

Relationship between GPA and other related methods
In this subsection, we briefly discuss the relationship of GPA with other related methods, such as LMM, conditional FDR and GSEA.
LMM is an effective tool for exploring genetic architecture of complex traits/diseases and it has been implemented in a popular software named GCTA [42]. Compared with LMM, GPA has the following distinct features: • LMM explores the genetic architecture underlying of complex trait/disease by estimating the gross phenotypic variance that can be explained by whole genome or a certain subset of SNPs. In contrast, GPA provides more "fine-grained" understanding by giving estimates of the local fasle discovery rate of each SNP, the proportion of SNPs that are associated with the phenotype ({π l }), the overall effect strength of the associated SNPs (α), and enrichment of a particular functional annotation ({q l }). • Application of LMM requires the availability of genotype data, while GPA only needs the summary statistics (p-values) as its input. Typically, the genotype data may not be accessible as easily as the summary statistics. For example, when researchers want to implement the integrative analysis for their own GWAS data at hand with related GWAS studies, it is much easier for them to obtain the summary statistics than the whole data sets of related GWAS studies. In this sense, GPA may greatly simplify the procedure of integrative analysis.
To our best knowledge, the conditional FDR approach is the first approach that statistically addresses the issue of pleiotropy between two GWAS, and GSEA is the most popular approach to evaluating the enrichment of gene sets. In fact, GPA provides a unified framework for systematically integrating both sources of information: pleiotropy and annotation. Rigorous statistical inference of pleiotropic effects and annotation enrichment has been established in this framework. As demonstrated in our simulation study (based on the widely accepted randomeffects model in GWAS), GPA has better performance of identifying disease-associated markers than the conditional FDR approach, and it shows a greater power of evaluating annotation enrichment than GSEA as well.

Conclusion
A large number of GWAS have been carried out in the last few years. As the cumulation of GWAS data, it is time to investigate systematic analysis of GWAS data sets to provide comprehensive understanding of the genetic architecture of complex traits/diseases, and provide new insights for functional genome studies. To achieve this goal, there have been fast growing interests in developing computational and statistical approaches to exploring genomic data in the post-GWAS era. In this paper, we present our approach, named GPA, which integrates information from multiple GWAS data sets and functional annotation data. Not only has GPA a better statistical power than related methods, it also provides interpretable model parameters which deepens our understanding of the genetic architecture of complex traits/diseases. We have successfully applied GPA to analyze GWAS data of five psychiatric disorders from PGC. The real data results demonstrate that GPA is able to identify pleiotropic effects among psychiatric disorders and detect enrichment of the CNS gene set. We have also applied GPA to analyze bladder cancer GWAS data with ENCODE data as annotation, where the significant enrichment of immune system and carcinoma pathways has been observed. In addition, GPA shows similar results of enrichment analysis to LMM, even without the genotype data. This suggests that GPA can serve as an effective tool for the integrative analysis of multiple GWAS data with functional annotation data, when genotype data are not available.
In this work, we mainly focused on studying the role of a single annotation vector. To deal with multiple annotation vectors, we assumed the conditional independence among them (see GPA model (7)). As a result, the current version of GPA may become unreliable in the presence of multiple highly correlated annotation vectors. Hence, when multiple sources of annotation data are available, we suggest incorporating them into GPA one at a time, as demonstrated in our analysis of bladder cancer GWAS data with ENCODE annotation. We realize that simultaneous analysis of multiple correlated annotation vectors may be more powerful. We will investigate this issue in our future work.

Web source
The R package of GPA is publicly available from GitHub http://dongjunchung.github.io/GPA/ and will be available in Bioconductor upon acceptance of the manuscript.
Based on this complete likelihood, the E-and M-steps in t-th iteration of the EM algorithm are obtained as follows.
. The M step remains the same when annotation data is not incorporated, except that we do not need to calculate q