Meta-analysis of set-based multiple phenotype association test based on GWAS summary statistics from different cohorts

Abstract

Genome-wide association studies (GWAS) have emerged as popular tools for identifying genetic variants that are associated with complex diseases. Standard analysis of a GWAS involves assessing the association between each variant and a disease. However, this approach suffers from limited reproducibility and difficulties in detecting multi-variant and pleiotropic effects. Although joint analysis of multiple phenotypes for GWAS can identify and interpret pleiotropic loci which are essential to understand pleiotropy in diseases and complex traits, most of the multiple phenotype association tests are designed for a single variant, resulting in much lower power, especially when their effect sizes are small and only their cumulative effect is associated with multiple phenotypes. To overcome these limitations, set-based multiple phenotype association tests have been developed to enhance statistical power and facilitate the identification and interpretation of pleiotropic regions. In this research, we propose a new method, named Meta-TOW-S, which conducts joint association tests between multiple phenotypes and a set of variants (such as variants in a gene) utilizing GWAS summary statistics from different cohorts. Our approach applies the set-based method that Tests for the effect of an Optimal Weighted combination of variants in a gene (TOW) and accounts for sample size differences across GWAS cohorts by employing the Cauchy combination method. Meta-TOW-S combines the advantages of set-based tests and multi-phenotype association tests, exhibiting computational efficiency and enabling analysis across multiple phenotypes while accommodating overlapping samples from different GWAS cohorts. To assess the performance of Meta-TOW-S, we develop a phenotype simulator package that encompasses a comprehensive simulation scheme capable of modeling multiple phenotypes and multiple variants, including noise structures and diverse correlation patterns among phenotypes. Simulation studies validate that Meta-TOW-S maintains a desirable Type I error rate. Further simulation under different scenarios shows that Meta-TOW-S can improve power compared with other existing meta-analysis methods. When applied to four psychiatric disorders summary data, Meta-TOW-S detects a greater number of significant genes.

1 Introduction

Genome-wide association study (GWAS) is typically employed to identify an individual genetic variant associated with a specific phenotype. However, in cases where causal variants have weak effects on the trait and it is challenging to detect these variants, set-based tests are employed to identify the joint effects of multiple variants (multi-variants) on a particular phenotype (Lin et al., 2022; Huang et al., 2011; Liu et al., 2010; Neale and Sham, 2004). Compared to single-variant approaches, set-based tests can help to reduce the number of genome-wide association tests and are effective when a causal variant is unobserved or multiple causal variants are present (Lin et al., 2022). Consequently, set-based tests have been applied in gene-set analyses of common variants and rare variants (Liu et al., 2014). For instance, both weighted and unweighted burden tests have a good performance when variants in a gene affect a phenotype in the same direction (Morris and Zeggini, 2010). The sequence kernel association test (SKAT) is a score-based variance-component test that accommodates variants in a gene with opposite effects on a phenotype (Wu et al., 2011). However, the performance of the approaches mentioned above depends on the weighting scheme, such as the MAF-based weighting scheme that up-weights the contribution of rare variants and down-weights that of common variants. This weighting scheme may lead to a loss in power when common variants in the region are associated with the phenotype (Dutta et al., 2019b; Dutta et al., 2019a). The optimal strategy for grouping or weighting genetic variants depends on the unknown genetic architecture of each phenotype and each variant. TOW for Testing the effect of an Optimal Weighted combination of variants in a set is a powerful method to increase power by assigning optimal weights to genetic variants (Sha et al., 2012).

Many complex phenotypes are influenced by multiple genetic variants, each with a small effect size. Set-based tests that consider a single phenotype might not capture the collective effects of these variants, leading to reduced power. Alternatively, cross-phenotype tests that aggregate associations in multiple phenotypes can substantially improve power over single phenotype-based methods (Dutta et al., 2019b; Dutta et al., 2019a). Meta-analysis of multiple phenotypes, using GWAS summary statistics, is a practical approach to increase power by increasing sample sizes and aggregating variants with small effect sizes to detect more significant pleiotropic genes (Cirulli et al., 2020; Panagiotou et al., 2013). Detecting the pleiotropic genes can provide insights into biological mechanisms influencing complex human phenotypes. The challenge in meta-analysis is that there is no uniformly most powerful (UMP) test. The power depends on signal directions and between-phenotype correlation. To boost analysis power, several methods have been proposed for GWAS multiple phenotype analysis. For instance, Fisher’s method of combining independent p-values has been extended to dependent univariate tests (Li et al., 2014). However, the p-value approximations of these tests are not accurate for small significance levels often required by GWASs. The minimum of the p-values (MinP) of multiple phenotypes has been proposed as a testing statistic (Conneely and Boehnke, 2007). And this test is powerful when a gene affects only a very small number of multiple phenotypes, but is less powerful in the presence of denser signals (Liu and Lin, 2018). The aggregated Cauchy association test (ACAT) is a flexible and computationally efficient p-value combination method that boosts power under various genetic architectures (Liu and Xie, 2019; Liu et al., 2019). Thus, in this article, we applied this approach to combine p-values of set-based tests from multiple cohort studies. Most importantly, the Cauchy combination is an extremely fast omnibus testing procedure that performs multiple testing adjustments analytically and applies to the combination of any tests (Li et al., 2023; Li et al., 2022; Li et al., 2020).

Meta-analysis is a promising approach to detecting a series of gene-phenotype associations that would have remained undetected by a participating cohort alone. In addition, a meta-analysis based on GWAS summary statistics simplifies data sharing, keeping sensitive individual data at the cohort level and sharing only non-sensitive summary data (Lin et al., 2022). For instance, Meta-MultiSKAT performs a variance component test and uses summary statistics to test for association between multiple continuous phenotypes and variants in a region. However, the p-value calculation of Meta-MultiSKAT relies on the normality assumption of the score vector and this assumption may be violated in the meta-analysis (Dutta et al., 2019b; Dutta et al., 2019a). MetaUSAT is a novel unified association test of multiple traits with only a single genetic variant, and the test statistic is dominated by a predefined parameter weight (Ray and Boehnke, 2018). In the cross meta-analysis, the overlapping subjects can induce a correlation between the summary statistics and inflate the false discovery rate of meta-analyses. Method FOLD is proposed to account for overlapping subjects at the summary statistics level using a split prior which categorizes subjects based on their contributions to the final statistic (Kim et al., 2017). However, it is designed for qualitative phenotypes and is difficult to obtain splitting prior if the numbers of cases and controls in any of the GWAS cohorts are missing. Here we propose a new approach, Meta-TOW-S, which conducts joint association tests between multiple phenotypes and genetic variants within a gene, utilizing GWAS summary statistics from diverse GWAS. Our approach applies set-based tests using an optimal weighted combination of variants and accounts for sample size differences across different GWAS by employing the Cauchy combination method. Meta-TOW-S combines the advantages of set-based tests and multi-phenotype modeling, exhibiting computational efficiency and enabling analysis across multiple phenotypes while accommodating overlapping samples from different cohorts.

To evaluate the performance of Meta-TOW-S, we need to mimic GWAS summary data from different cohorts. Thus, we also develop a phenotype simulator package that encompasses a comprehensive simulation scheme capable of modeling multiple phenotypes with multiple underlying genetic loci, intricate noise structures, and diverse correlation patterns among the phenotypes. The R package for the phenotype simulator is available on GitHub (https://github.com/Julia-lirong/PheGen). Furthermore, we evaluate the performance of our method using simulation studies and compare the power of our method with the power of three existing methods which integrate Burden (Morris and Zeggini, 2010), SKAT (Wu et al., 2011), and VEGAS (Liu et al., 2010) with Cauchy combination (Liu and Xie, 2019) to detect pleiotropic effects. Our simulation studies validate that Meta-TOW-S maintains a desirable Type I error rate and enhances power across various simulation scenarios compared with other existing meta-analysis methods. We also apply Meta-TOW-S to four psychiatric disorders summary data which are available from the Psychiatric Genomics Consortium (PGC) (Sullivan et al., 2018). The real data analyses demonstrate that Meta-TOW-S outperforms other comparison methods by detecting a greater number of significant genes.

2 Materials and methods

Consider phenotypes from GWAS cohorts with sample sizes that are subject to . For the cohort, suppose that the GWAS summary statistics are the Z-scores of genetic variants in a genomic region. We assume all cohorts share the same genetic variants in the specific region. is the number of overlapping subjects among all cohorts, where .

2.1 Meta-TOW-S

For the cohort, we suppose that is a matrix of genotypes in the interested genomic region (gene or pathway), and is a vector of phenotypes (either a quantitative or qualitative phenotype). In TOW (Sha et al., 2012), the generalized linear regression model with the fixed effect is used to link the phenotype and genotypes. The statistic model is defined as , where is the vector of weights for the genetic variants; , and is the weight assigned to the variant in the cohort. is the effect size of the weighted combination of genetic variants on the phenotype . Under the null hypothesis of no association between the variants in the region and the phenotype, we test . The score test statistic is given by , where , , , represents a vector containing all ones, and is a identity matrix. Under the null hypothesis, follows a multivariate normal distribution with mean vector and covariance matrix . The TOW method obtains the optimal weights by maximizing the score test statistic using the Cauchy-Schwartz inequality. Specifically, we have the form , with the equality when , where is any positive definite matrix, and is a constant. Based on Yan’s work (Yan, 2022), is a quadric term with the asymptotical distribution of weighted sum of Chi-squares. We assume is full rank and let . Then the optimal weight is obtained with the form . Using the optimal weights, the test statistic of TOW is given by . Under the null hypothesis, asymptotically follows a Chi-square distribution with a degree of freedom , that is (Yan, 2022; Sha et al., 2012).

Inspired by the idea proposed in this project (Svishcheva et al., 2019), Yan (2022) rewrite the test statistic using GWAS summary statistics. For the cohort, the Z score of the genetic variants in a region can be written as , where is an diagonal matrix of the square roots of genotypic variances of the variants. We assume that under the null hypothesis, follows the multivariate normal distribution , where is an matrix of correlations between the genotypes of these variant and . Then the test statistic of TOW based on individual data can be written as . This test statistic is called and asymptotically, where is the number of variants left after correlation pruning (Svishcheva et al., 2019). We can estimate using a reference sample of genotypes from the same population, such as 1,000 Genome phase 3 if individual genotype data are not available.

Consider phenotypes from different GWAS cohorts, we denote as the p-value of for the GWAS cohort, where . Then to detect the association between genetic variants in this region and multiple phenotypes, we define the Cauchy combination test statistic as , where the weight is defined as and has a standard Cauchy distribution under the null (Liu and Xie, 2019). Here we assigned more weight to a large GWAS cohort because a GWAS cohort with a large sample size carries more information than a smaller GWAS cohort (Zhu et al., 2015).

2.2 Comparison with other set-based association tests

Versatile set-based association study (VEGAS): For a specific region with genetic variants in the GWAS summary study, the test statistic of VEGAS is the sum of all squared Z-scores which is (Liu et al., 2010). Under the null hypothesis, asymptotically follows a mixture of chi-square distribution. To obtain the p-value of VEGAS, several methods have been proposed, such as numerical inversion of the characteristic function (Liu et al., 2009), Davies method (Davies, 1980), or Saddlepoint approximation (Kuonen, 1999).

SKAT and Burden: For the cohort, we consider the generalized linear model with the random effect, we have , where is a vector of effect size of the variants which is assumed to follow a normal distribution under the null hypothesis, where is the variance component of the phenotype explained by the variants. To test the association between the genetic variants in a region and the phenotype, Burden and SKAT test the null hypothesis against (Svishcheva et al., 2019). For the GWAS cohort, the burden test statistic using GWAS summary statistics can be written as , where is a diagonal matrix with , and is the minor allele frequency for the variant (Lee et al., 2013). Under the null hypothesis, follows a scaled Chi-square distribution with one degree of freedom. In SKAT, the test statistic based on GWAS summary statistic is defined as , where is a diagonal matrix of weights with . The beta distribution density function has pre-specified parameters , and . To obtain the p-value of Burden and SKAT, several methods have been proposed, such as numerical inversion of the characteristic function (Liu et al., 2009), Davies method (Davies, 1980), or Saddlepoint approximation (Kuonen, 1999).

Let be the p-value of the cohort study for based on the three comparison methods Vegas, Burden, and SKAT. We use the same strategy to detect the association between multiple phenotypes in different GWAS cohorts and genetic variants in a region by employing the Cauchy combination (Liu and Xie, 2019). And we designate these three methods as Meta-Vegas, Meta-Burden, and Meta-SKAT, respectively.

3 Simulation studies

3.1 Phenotype simulator

Suppose we have cohorts with sample sizes , respectively. For the cohort, the phenotype is modeled by a linear model , where are matrices of genotypes with columns standardized to mean zero and variance 1, with dimension . is a vector of standardized phenotypes with mean zero and variance 1. is the vector of genotypes effect sizes in the specific gene and is the vector of residuals representing environmental effects and non-additive genetic effects for the cohort. For each cohort, we assume that all genotype effect sizes are drawn with equal variance for all causal variants in a gene (Lee et al., 2014). Then for all cohorts, we suppose that the matrix of effect size has mean zero and covariance matrix , where is the genetic covariance on two phenotypes and on shared individuals, is the heritability explained by the variants in a region for the phenotype. The residuals has mean zero and covariance matrix , where is the covariance of non-genetic effects on the and phenotypes among shared individuals. Then, the phenotypic correlation for cohort and cohort among the overlapping samples is (Lemma). For these overlapping individuals, all cohorts share the same genotype data.

Next, to generate genotypes for individuals in a cohort, we employ the calibration coalescent model (COSI) to generate 10,000 haplotypes for a region of approximately 200 kbps, mimicking the LD structure found in individuals of European ancestry (He et al., 2017). We randomly select regions of 10 kbps in length, which encompass approximately 100 genetic variants, and utilize the simulated haplotypes to create genotypes for the variant sets. Among these genetic variants, we specifically designate 10% as causal variants, comprising 60% rare variants and 40% common variants, respectively. Subsequently, we generate the genetic component and non-genetic component based on the distribution described above. Specifically, we fix the phenotypic correlation between each pair of phenotypes. The phenotype correlation among the overlapping individuals is influenced by both genetic and non-genetic covariance, as proven in Lemma in the Supplementary Material. We allocate 80% of the phenotypic correlation to genetic covariance and 20% to non-genetic covariance. Finally, we generate K quantitative phenotypes in different cohorts using the additive model , .

3.2 Simulation

To set up a multi-cohort scenario, we generated individuals for multiple cohorts with different sample sizes but the same overlapping sample size for simplicity. To achieve a normally distributed input for the association test between gene and phenotype within each cohort, a rank-based inverse-normal transformation to the residuals of each phenotype was performed. In simulations, we access the performance of Meta-TOW-S with compared methods Meta-Vegas, Meta-Burden, and Meta-SKAT. We design different patterns of phenotypes by varying the correlation of phenotypes, sample sizes, and overlapping samples.

4 Results

4.1 Type I error rates

To evaluate the Type I error rates of Meta-TOW-S, we first simulate regions, each encompassing approximately 100 genetic variants, reflecting the LD structure observed in individuals of European ancestry. We then replicate times to generate the phenotypes under the null hypothesis of no genetic contribution to any of the three traits, that is . Then we simulate datasets to estimate the Type I error rates at nominal significance levels . We generate three phenotypes with different sample sizes. In the first situation, the sample sizes for the three phenotypes are equal, with a ratio of . In the second scenario, the sample sizes are unequal, with a ratio of . For simplicity, the number of overlapping individuals between any two phenotypes is kept constant at . The phenotypic correlation between any two phenotypes among the overlapping individuals is set to be either or . A correlation of 1 indicates simulation of mimicking one phenotype but form different cohorts

For Type I error evaluations, we use the significance levels . Table 1 summarizes the estimated Type I error rates of the four tests based on different settings. We can see from Table 1, that the Type I error rates of all tests are all within the estimated 95% confidence intervals in most situations indicating that the Type I error rates of the four tests are well controlled at the nominal significance levels.

4.2 Power comparisons

We compare the empirical power of Meta-TOW-S with Meta-Vegas, Meta-SKAT, and Meta-Burden. The power is defined as the proportion of test statistics with p-values less than the nominal significance level and we evaluate power at the nominal significance level after Bonferroni correction. For power comparisons, we generate phenotypes under the alternative hypothesis where the genetic effect is added correspondingly. Under each simulation setting, we generate 10,000 datasets to evaluate power at the nominal significance level . For the overlapping individuals, the phenotypic correlation between any pair of phenotypes across four scenarios: , and . A correlation of 1 indicates a simulation of mimicking one phenotype but from different cohorts. The primary simulations explore four distinct schemes. The first scheme compares the set-based association test for multiple phenotypes with that for a single phenotype. The second scheme evaluates the set-based association test for multiple phenotypes across varying numbers of traits, which are and . The third scheme tests the performance of a weighting scheme we designed within the Cauchy combination. The final scheme evaluates the set-based multiple phenotypes association test under different degrees of sample overlap. For simplicity, we set the heritability for each phenotype in the above scheme to be either or for . The simulation results are depicted in Figures 1–4 for a heritability of , while results for a heritability of are presented in the Supplementary Material.

FIGURE 1

FIGURE 2

FIGURE 3

FIGURE 4

We first compare the meta-analysis with single phenotype-based set-level tests. It shows that the integration of different cohorts to detect the genetic variants in a region associated with at least one phenotype could boost power (Figure 1; Supplementary Figure S1). Our proposed Meta-TOW-S has slightly better performance compared with Meta-Vegas, and Meta-SKAT, and is better performed than Meta-Burden since the Burden-based method has poor performance if the effects of causal variants are in different directions. We also found the increase in correlation between pairs of phenotypes will boost power as well, indicating that the integration of phenotypes will benefit to detect more pleiotropic genes.

Next, we vary the number of cohorts where each cohort has the same sample size and genetic heritability . Specifically, we create datasets with , , and cohorts and assess the effect of incorporating multiple cohorts in the meta-analysis for the identification of pleiotropic genes. As illustrated in Figure 2; Supplementary Figure S2, Meta-TOW-S outperforms or performs equivalently well as the other three methods. It also shows that the increase in the number of cohorts boosts the power of all methods.

In Meta-TOW-S, Meta-Vegas, Meta-SKAT, and Meta-Burden, more weight is assigned to a large study, and a small weight is assigned to a small study in the Cauchy combination. We compare those four methods with the unweighting scheme where the weight is the same among all studies in the Cauchy combination. Figure 3; Supplementary Figure S3 shows that the weighting scheme-based meta-analysis has a higher power compared with an unweighting scheme. Specifically, when the correlation between any pair of phenotypes is low, the weighting scheme meta-analysis exhibits significantly enhanced statistical power. Conversely, when the correlation between any pair of phenotypes is high, there is only a slight improvement in power. It indicates that the weighting scheme has a slight improvement when the genetic effect contributes to all three studies equivalently if the phenotype correlation is , which represents homogeneity across cohorts for the same phenotype.

Last, we consider the situation in which the sample sizes are the same in each cohort but different proportions of overlapping individuals are shared among all cohorts. As expected in Figure 4; Supplementary Figure S4, Meta-TOW-S outperforms the other three methods, and Meta-Vegas and Meta-SKAT have comparable performance. It shows that the power is further improved when there are larger overlapping samples between studies if phenotypes are highly correlated, which could be attributed to reduced heterogeneity across cohorts.

4.3 Real data analysis

We apply these four methods to GWAS summary statistics for four psychiatric disorders available from the Psychiatric Genomics Consortium (PGC) (Sullivan et al., 2018). These phenotypes are attention-deficit/hyperactivity disorder (ADHD), autism spectrum disorder (ASD), bipolar disorder (BD), and schizophrenia (SCZ) (Zhang et al., 2021). The sample sizes for these four traits range from to , with all individuals of European ancestry. We utilize the LD structure data from the 1,000 Genome Project Phase III European population as the reference in set-based multiple phenotype association test. The details of these four GWAS summary statistics are summarized in Supplementary Table S1. The analysis of these four methods for the four psychiatric disorders is summarized using the UpSet plot shown in Figure 5. We use the gene-based GWAS significance level in the analysis. For Meta-Burden and Meta-SKAT, the weight in the gene-based tests are defined in relation to the minor allele frequency (MAF) of genetic varaints. Specifically, we use the default weight for Meta-Burden and for Meta-SKAT, where the beta distribution density function has pre-specified parameters , , and of the variant. As a result, there are 557 genes detected by Meta-TOW-S, 517 genes detected by Meta-Vegas, 83 genes detected by Meta-SKAT, and 62 genes detected by Meta-Burden. These results indicate that Meta-TOW-S identifies more significant genes compared to Meta-Vegas, Meta-SKAT, and Meta-Burden.

FIGURE 5

4.4 Enrichment analysis

We implement the Gene Set Enrichment Analysis (GSEA) to analyze the significant genes identified by Meta-TOW-S that are enriched toward the top list of genes that are associated with specific biological pathways, processes, functions, or diseases (Subramanian et al., 2005). Gene Ontology (GO) is a community-based bioinformatics resource that employs ontologies to represent biological knowledge and describes information about gene and gene product information (Peng et al., 2017). Go is widely used to infer functional information for gene products, such as gene function enrichment, protein function prediction, and disease association analysis. And Go contains three categories: cellular component (CC), molecular function (MF: the biological function of the gene), and biological process (BP: pathways or larger processes that multiple gene products are involved in). KEGG (Kyoto Encyclopedia of Genes and Genomes) is a database resource for understanding high-level functions and utilities of the biological system. KEGG is used to search for the pathways associated with the identified genes. Detection of KEGG pathway database over-representation against a universal Homo Sapien background is assessed by hypergeometric tests (Solomon et al., 2022; Kanehisa et al., 2016). DisGeNET is an integrative and comprehensive resource of gene-disease associations from several public data sources and the literature (Piñero et al., 2015; Yu et al., 2014). It contains gene-disease associations and variant-gene-disease associations. The disease enrichment analysis is used to assess whether the genes associated with multiple phenotypes in meta-analysis are overrepresented in specific gene sets. A Bonferroni corrected cutoff of 0.05 was used for the significance of the pathway and disease. For genes annotation, we employed the org.Hs.e.g.,.db package in R. This package offers an extensive set of annotations for the human genome, including mappings between different gene identifiers and detailed genomic features.

We use the 557 significant genes identified by Meta-TOW-S, which are associated with four psychiatric disorders, for the enrichment analysis. We perform enrichment analyses at two different levels: pathways and diseases (Figures 6 and 7; Supplementary Figure S5). The top 20 enriched KEGG pathways are summarised in Figure 6 with results sorted from lowest to highest p-values. Of the 280 KEGG pathways tested, 37 were statistically significant after adjusting for multiple testing. The top KEGG pathways predominantly belong to groups associated with Human T-cell Leukemia virus 1 infection, Epstein-Barr virus infection, Antigen processing and presentation, etc. In the disease enrichment analysis, of the 5,496 disease tests, 455 had Bonferroni-corrected enrichment p-values lower than 0.05. The implicated genes are involved in Child Development Disorders Pervasive, Myasthenia Gravis, Vitiligo, Sarcoidosis, etc. (Figure 7).

FIGURE 6

FIGURE 7

5 Discussion

Much research suggests that many genes are associated with multiple correlated or even distinct phenotypes, and such associations have been termed cross-phenotype associations, which is relevant to pleiotropy in complex phenotypes (Li and Zhu, 2017). We propose a new method, Meta-TOW-S, which integrates association evidence of multiple phenotypes from study-specific GWAS summary statistics and thus detects the significant pleiotropic genes. This method Meta-TOW-S is based on the set-based test which uses the weights that maximize the score test statistic to increase power. To combine the test statistics from multiple GWAS cohorts, Meta-TOW-S uses the Cauchy combination by assigning more weight to a large GWAS to account for more biological information. In addition, Meta-TOW-S enables analysis across multiple phenotypes by accommodating overlapping samples from different cohorts.

To mimic real multiple comprehensive study-specific phenotypes, we develop a phenotype simulator that encompasses a simulation scheme capable of modeling multiple phenotypes with multiple underlying genetic loci, intricating noise structures, and different correlation patterns among the phenotypes with overlapped samples across different cohort studies. Our simulations show that the Type I error rates of Meta-TOW-S are well maintained under different conditions of phenotype correlation structures and overlapping samples and are more powerful than the other three comparison methods under most scenarios. We also find that the power of meta-analysis is significantly increased compared to the single phenotype set-based tests. A higher phenotype correlation, larger overlapping samples across multiple cohort studies, and more cohorts can increase power as well. We apply Meta-TOW-S to the summary statistics of four psychiatric disorders provided by the Psychiatric Genomics Consortium (PGC): attention-deficit/hyperactivity disorder (ADHD), autism spectrum disorder (ASD), bipolar disorder (BD), and schizophrenia (SCZ). As a result, 557 significant cross-phenotype associations are identified by Meta-TOW-S which is more than the number of genes identified by the other three methods. In the enrichment analysis, the gene sets of Child Development Disorders Pervasive are more enriched for genes associated with these four psychiatric disorders. In the KEGG pathway analysis, the significant genes identified by Meta-TOW-S showed notable enrichment in immune-related pathways rather than neurological processes. However, it is important to recognize that there is substantial evidence linking immune system involvement to neurological disorders. For instance, research has indicated that developmental disorders, such as Autism Spectrum Disorders (ASDs), can involve significant immune activity, including neuroinflammation, which plays a crucial role in the pathophysiology of these conditions (Vargas et al., 2005; Ashwood and Van de Water, 2004). Additionally, studies have demonstrated that elevated levels of regulatory T cells are associated with an increased risk of Attention-Deficit/Hyperactivity Disorder (ADHD) (Çetin et al., 2022), suggesting that immune dysregulation may contribute to the manifestation of neurological symptoms in certain contexts.

Meta-TOW-S has multiple advantages for identifying cross-phenotype associations. First of all, Meta-TOW-S can integrate information from multiple cohort studies to increase power and has the potential to detect more pleiotropic genes. Secondly, the test statistic has a standard Cauchy distribution under the null hypothesis and greatly reduces the computing time. Thirdly, this method is based on GWAS summary statistics from different cohort studies and GWAS summary statistics are more accessible than individual-level phenotype and genotype data The last point is that Meta-TOW-S gives more weight to the study with a larger sample size. Meanwhile, the developed phenotypes simulator can mimic complex structures in a meta-analysis which can be applied in other cross-phenotype analyses. However, Meta-TOW-S leverages information from correlated phenotypes to enhance its statistical power. Simulation results demonstrate that Meta-TOW-S surpasses other methods in power performance when phenotypic correlations are strong. However, when the correlation between phenotypes is weak, Meta-TOW-S may not outperform other methods, as it does not rely heavily on borrowing information from other phenotypes in such scenarios.

Currently, the framework of Meta-TOW-S needs to estimate the correlation matrix from a reference panel due to the unavailability of individual genotype data. However, the choice of the reference panel may influence the performance of Meta-TOW-S. The second challenge is that we use a correlation pruning procedure to ensure that the correlation matrix is fully ranked, which may drop some correlated variants in a gene.

In summary, the Meta-TOW-S method is a very useful method for detecting gene associations of multiple phenotypes from different GWAS cohorts. Meta-TOW-S has robust power and can handle different scenarios such as diverse phenotype correlation, and intricating cohort studies. The computational efficiency of Meta-TOW-S can also improve genetic discovery for hundreds of phenotypes across multiple GWAS cohorts in compliance with data privacy.

References

• 1

  AshwoodP.Van de WaterJ. (2004). Is autism an autoimmune disease?Autoimmun. Rev.3, 557–562. 10.1016/j.autrev.2004.07.036

  • CrossRef
  • Google Scholar

• 2

  ÇetinF. H.UçaryıLMAZH.UçarH. N.ArtaçH.GülerH. A.DuranS. A.et al (2022). Regulatory T cells in children with attention deficit hyperactivity disorder: a case-control study. J. Neuroimmunol.367, 577848. 10.1016/j.jneuroim.2022.577848

  • CrossRef
  • Google Scholar

• 3

  CirulliE. T.WhiteS.ReadR. W.ElhananG.MetcalfW. J.TanudjajaF.et al (2020). Genome-wide rare variant analysis for thousands of phenotypes in over 70,000 exomes from two cohorts. Nat. Commun.11, 542. 10.1038/s41467-020-14288-y

  • CrossRef
  • Google Scholar

• 4

  ConneelyK. N.BoehnkeM. (2007). So many correlated tests, so little time! Rapid adjustment of P values for multiple correlated tests. Am. J. Hum. Genet.81, 1158–1168. 10.1086/522036

  • CrossRef
  • Google Scholar

• 5

  DaviesR. B. (1980). Algorithm AS 155: the distribution of a linear combination of χ 2 random variables. J. R. Stat. Soc. Ser. C Appl. Statistics29, 323–333. 10.2307/2346911

  • CrossRef
  • Google Scholar

• 6

  DemontisD.WaltersR. K.MartinJ.MattheisenM.AlsT. D.AgerboE.et al (2019). Discovery of the first genome-wide significant risk loci for attention deficit/hyperactivity disorder. Nat. Genet.51, 63–75. 10.1038/s41588-018-0269-7

  • CrossRef
  • Google Scholar

• 7

  DuttaD.Gagliano TaliunS. A.WeinstockJ. S.ZawistowskiM.SidoreC.FritscheL. G.et al (2019a). Meta‐MultiSKAT: multiple phenotype meta‐analysis for region‐based association test. Genet. Epidemiol.43, 800–814. 10.1002/gepi.22248

  • CrossRef
  • Google Scholar

• 8

  DuttaD.ScottL.BoehnkeM.LeeS. (2019b). Multi‐SKAT: general framework to test for rare‐variant association with multiple phenotypes. Genet. Epidemiol.43, 4–23. 10.1002/gepi.22156

  • CrossRef
  • Google Scholar

• 9

  GroveJ.RipkeS.AlsT. D.MattheisenM.WaltersR. K.WonH.et al (2019). Identification of common genetic risk variants for autism spectrum disorder. Nat. Genet.51, 431–444. 10.1038/s41588-019-0344-8

  • CrossRef
  • Google Scholar

• 10

  HeZ.XuB.LeeS.Ionita-LazaI. (2017). Unified sequence-based association tests allowing for multiple functional annotations and meta-analysis of noncoding variation in metabochip data. Am. J. Hum. Genet.101, 340–352. 10.1016/j.ajhg.2017.07.011

  • CrossRef
  • Google Scholar

• 11

  HuangH.ChandaP.AlonsoA.BaderJ. S.ArkingD. E. (2011). Gene-based tests of association. PLoS Genet.7, e1002177. 10.1371/journal.pgen.1002177

  • CrossRef
  • Google Scholar

• 12

  KanehisaM.SatoY.KawashimaM.FurumichiM.TanabeM. (2016). KEGG as a reference resource for gene and protein annotation. Nucleic acids Res.44, D457–D462. 10.1093/nar/gkv1070

  • CrossRef
  • Google Scholar

• 13

  KimE. E.LeeS.LeeC. H.OhH.SongK.HanB. (2017). FOLD: a method to optimize power in meta-analysis of genetic association studies with overlapping subjects. Bioinformatics33, 3947–3954. 10.1093/bioinformatics/btx463

  • CrossRef
  • Google Scholar

• 14

  KuonenD. (1999). Miscellanea. Saddlepoint approximations for distributions of quadratic forms in normal variables. Biometrika86, 929–935. 10.1093/biomet/86.4.929

  • CrossRef
  • Google Scholar

• 15

  LeeS.AbecasisG. R.BoehnkeM.LinX. (2014). Rare-variant association analysis: study designs and statistical tests. Am. J. Hum. Genet.95, 5–23. 10.1016/j.ajhg.2014.06.009

  • CrossRef
  • Google Scholar

• 16

  LeeS.TeslovichT. M.BoehnkeM.LinX. (2013). General framework for meta-analysis of rare variants in sequencing association studies. Am. J. Hum. Genet.93, 42–53. 10.1016/j.ajhg.2013.05.010

  • CrossRef
  • Google Scholar

• 17

  LinA.ShadrinA.Van der MeerD.HindleyG.ChengW.SønderbyI. E.et al (2022). Efficient meta-analysis of multivariate genome-wide association studies with Meta-MOSTest. bioRxiv.

  • Google Scholar

• 18

  LiQ.HuJ.DingJ.ZhengG. (2014). Fisher's method of combining dependent statistics using generalizations of the gamma distribution with applications to genetic pleiotropic associations. Biostatistics15, 284–295. 10.1093/biostatistics/kxt045

  • CrossRef
  • Google Scholar

• 19

  LiuD. J.PelosoG. M.ZhanX.HolmenO. L.ZawistowskiM.FengS.et al (2014). Meta-analysis of gene-level tests for rare variant association. Nat. Genet.46, 200–204. 10.1038/ng.2852

  • CrossRef
  • Google Scholar

• 20

  LiuJ. Z.McraeA. F.NyholtD. R.MedlandS. E.WrayN. R.BrownK. M.et al (2010). A versatile gene-based test for genome-wide association studies. Am. J. Hum. Genet.87, 139–145. 10.1016/j.ajhg.2010.06.009

  • CrossRef
  • Google Scholar

• 21

  LiuH.TangY.ZhangH. H. (2009). A new chi-square approximation to the distribution of non-negative definite quadratic forms in non-central normal variables. Comput. Statistics & Data Analysis53, 853–856. 10.1016/j.csda.2008.11.025

  • CrossRef
  • Google Scholar

• 22

  LiuY.ChenS.LiZ.MorrisonA. C.BoerwinkleE.LinX. (2019). ACAT: a fast and powerful p value combination method for rare-variant analysis in sequencing studies. Am. J. Hum. Genet.104, 410–421. 10.1016/j.ajhg.2019.01.002

  • CrossRef
  • Google Scholar

• 23

  LiuY.XieJ. (2019). Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures. J. Am. Stat. Assoc.115, 393–402. 10.1080/01621459.2018.1554485

  • CrossRef
  • Google Scholar

• 24

  LiuZ.LinX. (2018). Multiple phenotype association tests using summary statistics in genome‐wide association studies. Biometrics74, 165–175. 10.1111/biom.12735

  • CrossRef
  • Google Scholar

• 25

  LiX.LiZ.ZhouH.GaynorS. M.LiuY.ChenH.et al (2020). Dynamic incorporation of multiple in silico functional annotations empowers rare variant association analysis of large whole-genome sequencing studies at scale. Nat. Genet.52, 969–983. 10.1038/s41588-020-0676-4

  • CrossRef
  • Google Scholar

• 26

  LiX.QuickC.ZhouH.GaynorS. M.LiuY.ChenH.et al (2023). Powerful, scalable and resource-efficient meta-analysis of rare variant associations in large whole genome sequencing studies. Nat. Genet.55, 154–164. 10.1038/s41588-022-01225-6

  • CrossRef
  • Google Scholar

• 27

  LiX.ZhuX. (2017). Cross-phenotype association analysis using summary statistics from GWAS. Stat. Hum. Genet. Methods Protoc.1666, 455–467. 10.1007/978-1-4939-7274-6_22

  • CrossRef
  • Google Scholar

• 28

  LiZ.LiX.ZhouH.GaynorS. M.SelvarajM. S.ArapoglouT.et al (2022). A framework for detecting noncoding rare-variant associations of large-scale whole-genome sequencing studies. Nat. methods19, 1599–1611. 10.1038/s41592-022-01640-x

  • CrossRef
  • Google Scholar

• 29

  MorrisA. P.ZegginiE. (2010). An evaluation of statistical approaches to rare variant analysis in genetic association studies. Genet. Epidemiol.34, 188–193. 10.1002/gepi.20450

  • CrossRef
  • Google Scholar

• 30

  NealeB. M.ShamP. C. (2004). The future of association studies: gene-based analysis and replication. Am. J. Hum. Genet.75, 353–362. 10.1086/423901

  • CrossRef
  • Google Scholar

• 31

  PanagiotouO. A.WillerC. J.HirschhornJ. N.IoannidisJ. P. (2013). The power of meta-analysis in genome-wide association studies. Annu. Rev. genomics Hum. Genet.14, 441–465. 10.1146/annurev-genom-091212-153520

  • CrossRef
  • Google Scholar

• 32

  PardiñasA. F.HolmansP.PocklingtonA. J.Escott-PriceV.RipkeS.CarreraN.et al (2018). Common schizophrenia alleles are enriched in mutation-intolerant genes and in regions under strong background selection. Nat. Genet.50, 381–389. 10.1038/s41588-018-0059-2

  • CrossRef
  • Google Scholar

• 33

  PengJ.WangH.LuJ.HuiW.WangY.ShangX. (2017). Identifying term relations cross different gene ontology categories. BMC Bioinforma.18, 573–574. 10.1186/s12859-017-1959-3

  • CrossRef
  • Google Scholar

• 34

  PiñeroJ.Queralt-RosinachN.BravoA.Deu-PonsJ.Bauer-MehrenA.BaronM.et al (2015). DisGeNET: a discovery platform for the dynamical exploration of human diseases and their genes. Database2015, bav028. 10.1093/database/bav028

  • CrossRef
  • Google Scholar

• 35

  RayD.BoehnkeM. (2018). Methods for meta‐analysis of multiple traits using GWAS summary statistics. Genet. Epidemiol.42, 134–145. 10.1002/gepi.22105

  • CrossRef
  • Google Scholar

• 36

  ShaQ.WangX.WangX.ZhangS. (2012). Detecting association of rare and common variants by testing an optimally weighted combination of variants. Genet. Epidemiol.36, 561–571. 10.1002/gepi.21649

  • CrossRef
  • Google Scholar

• 37

  SolomonO.HuenK.YousefiP.KüpersL. K.GonzálezJ. R.SudermanM.et al (2022). Meta-analysis of epigenome-wide association studies in newborns and children show widespread sex differences in blood DNA methylation. Mutat. Research/Reviews Mutat. Res.789, 108415. 10.1016/j.mrrev.2022.108415

  • CrossRef
  • Google Scholar

• 38

  StahlE. A.BreenG.ForstnerA. J.McquillinA.RipkeS.TrubetskoyV.et al (2019). Genome-wide association study identifies 30 loci associated with bipolar disorder. Nat. Genet.51, 793–803. 10.1038/s41588-019-0397-8

  • CrossRef
  • Google Scholar

• 39

  SubramanianA.TamayoP.MoothaV. K.MukherjeeS.EbertB. L.GilletteM. A.et al2005. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles, Proc. Natl. Acad. Sci. U. S. A.,102, 15545–15550. 10.1073/pnas.0506580102

  • CrossRef
  • Google Scholar

• 40

  SullivanP. F.AgrawalA.BulikC. M.AndreassenO. A.BørglumA. D.BreenG.et al (2018). Psychiatric genomics: an update and an agenda. Am. J. Psychiatry175, 15–27. 10.1176/appi.ajp.2017.17030283

  • CrossRef
  • Google Scholar

• 41

  SvishchevaG. R.BelonogovaN. M.ZorkoltsevaI. V.KirichenkoA. V.AxenovichT. I. (2019). Gene-based association tests using GWAS summary statistics. Bioinformatics35, 3701–3708. 10.1093/bioinformatics/btz172

  • CrossRef
  • Google Scholar

• 42

  VargasD. L.NascimbeneC.KrishnanC.ZimmermanA. W.PardoC. A. (2005). Neuroglial activation and neuroinflammation in the brain of patients with autism. Ann. Neurology Official J. Am. Neurological Assoc. Child Neurology Soc.57, 67–81. 10.1002/ana.20315

  • CrossRef
  • Google Scholar

• 43

  WuM. C.LeeS.CaiT.LiY.BoehnkeM.LinX. (2011). Rare-variant association testing for sequencing data with the sequence kernel association test. Am. J. Hum. Genet.89, 82–93. 10.1016/j.ajhg.2011.05.029

  • CrossRef
  • Google Scholar

• 44

  YanS. (2022). Statistical methods for controlling population stratification and gene-based association studies. Michigan Technological University.

  • Google Scholar

• 45

  YuG.WangL.-G.YanG.-R.HeQ.-Y. (2014). DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics31, 608–609. 10.1093/bioinformatics/btu684

  • CrossRef
  • Google Scholar

• 46

  ZhangY.LuQ.YeY.HuangK.LiuW.WuY.et al (2021). SUPERGNOVA: local genetic correlation analysis reveals heterogeneous etiologic sharing of complex traits. Genome Biol.22, 262–330. 10.1186/s13059-021-02478-w

  • CrossRef
  • Google Scholar

• 47

  ZhuX.FengT.TayoB. O.LiangJ.YoungJ. H.FranceschiniN.et al (2015). Meta-analysis of correlated traits via summary statistics from GWASs with an application in hypertension. Am. J. Hum. Genet.96, 21–36. 10.1016/j.ajhg.2014.11.011

  • CrossRef
  • Google Scholar