Here, we discuss the proposed meta-analysis framework, ancMETA. It performs meta-analysis at two different levels by aggregating multiple independent population GWAS summary statistics datasets. It uses an integrative analysis through Bayesian posterior probability and combines the results into known biological protein–protein network datasets. Lastly, ancMETA performs the meta-analysis at the subnetwork level and identifies the most significant subnetworks to understand the biological pathways (Figure 1).

We describe six different steps of the proposed meta-analysis framework as follows:

Step 1: Collection of N-independent studies.

This step requires N ≥ 2 independent GWAS summary statistics datasets, from the same phenotype or trait. GWAS summary statistics are defined here as per-genetic locus effect sizes (log odds ratios) together with their standard errors, p-values, or z-scores for the affected–unaffected traits de Leeuw et al. (2015); Il-Youp and Wei (2015); Huang et al., 2016); Wang et al. (2017).

Step 2: Mapping SNPs to the associated genes.

This is an intermediate step, and it is similar to our previous approach Chimusa et al. (2015), where all the SNPs are mapped to their related genes. It is common practice to assign SNPs to the genes based on a distance cutoff, and the previous studies use a variety of cutoffs, such as distance from 2 to 500 Kbps Conti et al. (2009); Brodie et al. (2016); Chimusa et al. (2015); Wang et al. (2010, 2007). ancMETA allows users to specify the distance cutoffs of assigning SNPs down-/upstream to a specific gene (see the Supplementary Material). At this step, the combined statistical outcomes (i.e., the effect size and the standard error) are computed, as illustrated in Figure 1.

We assume that the ith study (i = 1, 2,…, N) has J _i genes G j i ( j = 1,2 , … , J i ) , each having a set of specific SNPs, where N is the total number of studies. Let S = ∪ i , j S j i denote the set of all SNPs, where S j i is the set of SNPs at the jth gene of the study i. For an SNP s ∈ S j i , let β j i ( s ) , ζ j i ( s ) , and p j i ( s ) represent the effect size, the standard error, and the p-value, respectively. The proposed meta-analysis framework assumes that the SNPs related to one gene are correlated. This assumption is due to the fact that each gene contains a large number of SNPs mostly under LD; therefore, the magnitude of the statistical outcomes for each SNP within a gene from each study is considered to be approximately the same. Thus, the effect size and the p-value will follow a normal distribution. Let μ j i be the unknown true effect size at the jth gene G j i of study i (i = 1, 2,…, N; j = 1, 2,…, J _i ). Thus, a fixed-effects model can be applied to estimate the combined effect size at the gene level from each study i = 1, 2, … N. We construct a linear function of parameters μ j i as follows: x = y β , μ = ∑ s ∈ S j i β j i s μ j i .

Let V j i ( s ) be the variance. Then, the estimation of the combined effect size is given as follows: p x | B j i , μ j i , V j i = N x | B j i μ j i   , V j i ∼ ∏ s ∈ S j i 1 2 V j i s π exp − β j i s − μ j i 2 2 V j i s . (1)

Note Eq. 1 is the likelihood when μ ≠ 0. Although the solution, which corresponds to the maximum true effect size μ j i at gene G j i , can be approximated from sampling approaches, Eq. 1 can be solved analytically by differentiating the log likelihood and setting it to 0, that is, ∂ ⁡ log p x | B j i , μ j i , V j i ∂ μ j i = 0

or ∑ s ∈ S j i β j i s − μ j i V j i s = 0 . (2)

Eq. 2 can be written as ∑ s ∈ S j i V j i s − 1 β j i s − ∑ s ∈ S j i V j i s − 1 μ j i = 0 . Hence, the estimated effect size of gene G j i is given as μ ̂ j i = ∑ s ∈ S j i V j i s − 1 β j i s ∑ s ∈ S j i V j i s − 1 . Therefore, its standard error is derived as s d μ ̂ j i = ∑ s ∈ S j i V j i s − 1 .

Step 3: Meta-analysis at the gene level.

Let β j i and ζ j i be the combined statistical outcomes (i.e., effect size and standard error, respectively). Since the effect size may be different across the studies, the heterogeneity between the studies may be high. A parameter σ j 2 , which accounts for between-study variability, is introduced. Because the maximum-likelihood method underestimates the variance when the number of studies is small DerSimonian and Laird (1986), we use the unbiased method of moments Borenstein et al. (2010) to compute the variance components as follows: σ j 2 = 1 N ∑ i = 1 N β j i − H j 2 ,

where H j = 1 N ∑ i = 1 N β j i . It is assumed that the true effect size of the gene G _j is μ j + σ j 2 , where μ _j is an unknown parameter. The same calculations as in step 2 are applied to estimate the effect size of the gene G _j across the studies, which is calculated as follows: μ ̂ j = ∑ i = 1 N σ j 2 + V j i − 1 β j i ∑ i = 1 N σ j 2 + V j i − 1 ,

and its standard error is given as follows: s d μ ̂ j = ∑ i = 1 N σ j 2 + V j i − 1 .

If the precision ρ ⁱ of the study ith for i = 1, 2,…, N is equal to 1 ( μ i ) 2 , then the heterogeneity statistic, using Cochrane’s test (Q) is given as follows: Q = 100 × ∑ i N ρ i μ ̂ j − θ 2 , (3) where θ = Σ i N ρ i μ ̂ j Σ i N ρ i is a weighted estimator, and Q ∼ χ N − 1 2 under the null hypothesis that the between-study variance τ ² = 0. If we consider the expectation of Eq. 3, the between-study variance is given as τ 2 = Q − N − 1 Σ i N ρ i − Σ i N ρ i 2 Σ i N ρ i .

Since the test statistics in Eq. 3 has N − 1 degree of freedom (df), then for Q < df (which means τ ² < 0), the maximum between 0 and τ ² is considered, so that τ ² is non-negative. Recalling that in step 2, we pointed out that SNPs within one gene are correlated. Therefore, to perform meta-analysis, a fixed-effects model is used Borenstein et al. (2010).

Step 4: Mapping genes to biological networks.

For this step, all the genes and their related statistical outcomes and gene–gene combined LD are mapped into a PPI network database (Supplementary Material), which contains information on interactions among the genes. Therefore, the genes are considered as weighted nodes with their related statistical outcomes (i.e., effect size and standard error), and the edges represent the interactions between the genes (nodes), which are weighted with the combined LD, for each pair-wise genes (nodes) that have a link. More details can be found in Supplementary Materials. The details of the combined LD can be found in Supplementary Text S1.

Step 5: Meta-analysis at the subnetwork level. This step estimates the effect size and standard error of each subnetwork from . Supplementary Text S1 provides more details of different concepts used in . The procedure is the same as that of step 3, where it is assumed that there is a heterogeneity between the genes in each subnetwork. Let ρ j 2 be the between-gene heterogeneity, β j i be the observed effect size, and ζ j i be its standard error. If μ is the unknown true effect size at the subnetwork level with K genes, the estimated effect size is given by p ( x | β j i , μ j i , V j i + ρ j 2 ) ∼ ∏ G j i ∈ K 1 2 V j i + ρ j 2 π e − β j i − μ j i 2 2 V j i + ρ j 2 . (4)

Since Eq. 4 is considered to be the likelihood when μ ≠ 0, maximizing the log likelihood by solving ∂ ⁡ log p t | β j i , μ , V j i + ρ j 2 ∂ μ j i = 0

yields μ ̂ j i = ∑ G j i ∈ K ρ j 2 + V j i − 1 β j ∑ G j i ∈ K ρ j 2 + V j i − 1 , (5) and its standard error estimated effect size is given as s d μ ̂ j i = ∑ G j i ∈ K ρ j 2 + V j i − 1 .

Step 6: Computing the overall statistical significance.

From step 5, we know that β j i and ζ j i represent the observed effect size and the standard error, respectively, where μ ̂ j i as the unknown effect was estimated using Equation 5. For this step, let us consider the Bayes’ theorem, which includes all the possible effect outcomes for N studies, and let us apply the symmetrical property of the Gaussian distribution. Therefore, the posterior probability is approximately given as e ̂ j i = π ∫ R N β j i | μ ̂ j i , ζ j i p μ ̂ j i d μ ̂ j i 1 − π N β j i | 0 , ζ j i + π ∫ R N β j i | μ ̂ j i , ζ j i p μ ̂ j i d μ ̂ j i = π N β j i | X ̂ j i , ζ j i + V ̂ j i 1 − π N μ j i | 0 , ζ j i + π N μ j i | X ̂ j i , ζ j i + V ̂ j i . (6)

However, a test statistic based on the weighted sum of chi-square can be used, under the null hypothesis that there is no association. Thus, from the fact that e j i contains the information from all studies, it can be used as a prior weight in the weighted sum of the chi-square statistics as follows: χ ̂ j = ∑ i = 1 N e ̂ j i e ̂ j i μ ̂ j 2 ∑ i = 1 N ζ j i + τ 2 , τ 2 ≥ 0 , where τ ² is the between-study variance. Recalling that Equation 6 is not normally distributed, the p-values can be evaluated through a sampling method developed by Han and Eskin (2011). puts together all the steps of the proposed framework.

The U.S. residents of northern and western European ancestry (CEU), Yoruba (YRI), and Mexican (MEX) populations from the HapMap3 Project were used to generate independent case–control studies. Details of these populations can be found in Supplementary Table S1. We independently performed a population growth model on each dataset mentioned previously (Supplementary Text S2). From the resulting expanded datasets, we generated three independent case–control datasets based on chromosomes 1 and 22 using HapGen2 Su et al. (2011). Here, we randomly selected 3 SNPs on chromosome 1 and 3 other SNPs on chromosome 22 to be simulated as causal disease SNPs with differing effect sizes under HapGen2 Su et al. (2011). The simulated disease effect size parameters of those SNPs are summarized in Supplementary Table S2. These parameters are chosen to fit small effect size in some studies and strong in others. We simulated 1,000, 3,000, and 950 cases and 1,000, 3,000, and 1,000 controls from the haplotype (combination of multiple SNPs) data of CEU, YRI, and MEX, respectively. We considered the simulated MEX GWAS dataset as our primary study. The simulation details can be found in Supplementary Text S2.

We conducted a GWAS on each dataset using EMMAX Kang et al. (2010). As expected and according to our simulation parameters (Supplementary Table S2), the GWAS results indicate significant, moderate, and weak signals of association in CEU, YRI, and MEX, respectively, as per simulation (Supplementary Table S3 and Supplementary Figure S1). We used the resulting GWAS summary statistics as input in our proposed model, implemented in “ancMETA,” to perform the gene- and subnetwork-specific meta-analysis (Supplementary Text S3).

At the gene level, the results in Table 1 and Supplementary Figure S2 show significant genes, including CBX7, LD0C1L, and ASTN1 associated with our simulated causal SNPs (Supplementary Table S2). This result indicates further increase in the effect size across these 3 simulated case–control studies, compared to single-SNP tests in GWAS based on EMMAX (Supplementary Table S3). Interestingly, at a subnetwork level, we observed a significant convergence of effect sizes of the simulated disease variants across all the studies (Figures 2A,B), particularly to CBX7 (Table 1). This supports the fact that the true effect risk-associated variants may differ across populations at the polymorphism level, but the effect may tend to convert at a similar magnitude at gene and pathway levels or overlap in the same biological subnetworks/pathways when aggregating several small polymorphism effects.

Here, we used seven European bipolar disorder (BD) GWAS summary statistics obtained from the NIMH data repository (Supplementary Table S1). These datasets include Irish (IRI), Scottish (SCT), 3 European American (EUA), Norwegian (NOR), and British (BRB) individuals. Supplementary Table S4 provides GWAS summary results of each of these studies. Leveraging the European reference panel from 1000 Genomes Project data, we first conducted the imputation LD fine-mapping from the obtained GWAS summary statistics using ImpG-Summary Pasaniuc et al. (1977) to possibly unravel more LD-SNP association (Supplementary Figure S3). The resulting imputation GWAS LD fine-mapping summary statistics were used as input for ancMETA. We performed the meta-analysis at a gene level across these 7 studies, and seven genes were identified (AGT CACNA1C, ESR1, NCAN, BDNF, BCR, and GSK3B), of which AGT has a strong effect across the European populations including Irish, Scottish, Norwegian, and British in contrast to the European American (Supplementary Figures S4,5 and Supplementary Table S5).

Leveraging the recent version of the human PPI network (IntAct release 239) from the IntAct database Kerrien et al. (2012); Orchard et al. (2014), we observed a significant connected subnetwork where the estrogen receptor 1 (ESR1) gene is the main hub. ESR1 interacts with AGT and is connected with the other gene hubs with known BD-associated genes such as CACNA1C, PLCG2, NCAN, BCR, and BDNF. Supplementary Table S5 summarizes the association of the top subnetwork across these 7 studies. The identified subnetwork (Figures 3A,B) is significantly associated with the neurotrophin signaling (p = 4e ⁻¹⁴) and other interesting biological pathways such as myometrial relaxation and contraction (p = 3e ⁻⁰⁸), morphine addiction (p = 9.2e ⁻⁰⁸), DNA damage response (p = 3e ⁻⁰⁵), and alcoholism (p = 2e ⁻⁰⁵). In addition, the subnetwork (Figure 3C) is implicated in a positive regulation of transcription from RNA polymerase II promoter involved in neuronal differentiation (p = 7.7e ⁻²¹) and is also associated with autosomal dominant (p = 2.8e ⁻⁰⁷) inheritance types of diseases (Figure 3C).

To evaluate the determined interconnectivity among AGT, CACNA1C, ESR1, NCAN, BDNF, BCR, and GSK3B as well as with the other protein-coding genes in Figure 3, and their association with BD, we leverage GeneMania Franz et al. (2018) to reconstruct the network shown in Figure 3 based on physical and co-expressed interactions. The result from the reconstructed co-expression and physical interaction network of genes in Figure 3A confirms that the above 7 genes detected by ancMETA (Figure 3A) are physically and co-expressly interconnected Figure 4 and as well as with other protein-coding genes in Figure 3A. From the enrichment analysis based on Enrichr Kuleshov et al. (2016), the network in Figure 4 is significantly associated with neurotrophin signaling and other interesting biological pathways such as ErbB signaling pathways (Figure 4). Interestingly, recent studies have revealed that complex ErbB signaling networks regulate the assembly of neural circuitry, myelination, neurotransmission, and synaptic plasticity Mei and Nave (2014). Evidence indicates that there is an optimal level of ErbB signaling in the brain, and a deviation from it impairs brain functions. The ErbB signaling pathway may provide therapeutic targets for specific neuropsychiatric symptoms, and dysregulation in the ErbB signaling pathway may explain abnormalities of neural precursor migration in BD Mei and Nave (2014).