For the k t h cohort, we suppose that G k is a n k × M matrix of genotypes in the interested genomic region (gene or pathway), and y k is a n k × 1 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 g E y k | G k = β k 0 + G k w k β k , where w k is the vector of weights for the M genetic variants; w k = w k 1 , … , w k j , … , w k M T , and w k j is the weight assigned to the j t h variant in the k t h cohort. β k is the effect size of the weighted combination of genetic variants G k w k on the phenotype y k . Under the null hypothesis of no association between the variants in the region and the k t h phenotype, we test H 0 ∶ β k = 0 . The score test statistic is given by T k = n k w k T G k T P 0 y k y k T P 0 G k w k y k T P 0 y k w k T G k T P 0 G k w k = u k T w k w k T u k w k T Σ k w k , where P 0 = I n k − 1 n k 1 n k 1 n k T , u k = G k T P 0 y k / y k T P 0 y k / n k , Σ k = G k T P 0 G k , 1 n k represents a n k × 1 vector containing all ones, and I n k is a n k × n k identity matrix. Under the null hypothesis, u k follows a multivariate normal distribution with mean vector 0 and covariance matrix Σ k . The TOW method obtains the optimal weights w k * by maximizing the score test statistic T k using the Cauchy-Schwartz inequality. Specifically, we have the form T k = u k T w k w k T u k w k T Σ k w k = u k T E − 1 / 2 E 1 / 2 w k w k T E − 1 / 2 E 1 / 2 u k w k T Σ k w k ≤ E − 1 2 u k , E − 1 2 u k E 1 2 w k , E 1 2 w k w k T Σ k w k = w k T E w k w k T Σ k w k u k E − 1 u k = c u k E − 1 u k , with the equality when w k ∝ E − 1 u k , where E is any M × M positive definite matrix, and c = w k T E w k w k T Σ k w k is a constant. Based on Yan’s work (Yan, 2022), T k = c u k E − 1 u k is a quadric term with the asymptotical distribution of weighted sum of Chi-squares. We assume Σ k is full rank and let E = Σ k . Then the optimal weight is obtained with the form w k * = Σ k − 1 u k . Using the optimal weights, the test statistic of TOW is given by T k T O W = u k T ∑ k − 1 u k = n k y k T P 0 G k Σ k − 1 G k T P 0 y k y k T P 0 y k . Under the null hypothesis, T k T O W asymptotically follows a Chi-square distribution with a degree of freedom M , that is T k T O W ∼ χ M 2 (Yan, 2022; Sha et al., 2012).

Inspired by the idea proposed in this project (Svishcheva et al., 2019), Yan (2022) rewrite the test statistic T k T O W using GWAS summary statistics. For the k t h cohort, the Z score of the M genetic variants in a region can be written as S k = Z k 1 , … , Z k M T = V k − 1 G k T P 0 y k y k T P 0 y k / n k , where V k is an M × M diagonal matrix of the square roots of genotypic variances of the M variants. We assume that under the null hypothesis, S k follows the multivariate normal distribution N 0 , U , where U is an M × M matrix of correlations between the genotypes of these variant and U = V k − 1 G k T P 0 G k V − 1 n k . Then the test statistic of TOW based on individual data can be written as T k T O W = n k y k T P 0 G k Σ k − 1 G k T P 0 y k y k T P 0 y k = S k T U − 1 S k . This test statistic is called T T O W − S S k and T T O W − S S k ∼ χ M ′ 2 asymptotically, where M ′ is the number of variants left after correlation pruning (Svishcheva et al., 2019). We can estimate U 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 K phenotypes from different GWAS cohorts, we denote p k as the p-value of T T O W − S for the k t h GWAS cohort, where k = 1 , … , K . Then to detect the association between genetic variants in this region and multiple phenotypes, we define the Cauchy combination test statistic as T C = ∑ K k = 1 ν k ⁡ tan 0.5 − p k π , where the weight is defined as ν k = n k N and T C 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).

Versatile set-based association study (VEGAS): For a specific region with M genetic variants in the k t h GWAS summary study, the test statistic of VEGAS is the sum of all squared Z-scores which is T V e g a s S k = S k T S k (Liu et al., 2010). Under the null hypothesis, T V e g a s S k 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 k t h cohort, we consider the generalized linear model with the random effect, we have g E y k | G k = β k 0 + G k β k , where β k is a vector of effect size of the M variants which is assumed to follow a normal distribution N 0 , τ k 2 I under the null hypothesis, where τ k 2 is the variance component of the phenotype explained by the M variants. To test the association between the M genetic variants in a region and the phenotype, Burden and SKAT test the null hypothesis H 0 ∶ τ k 2 = 0 against H a ∶ τ k 2 > 0 (Svishcheva et al., 2019). For the k t h GWAS cohort, the burden test statistic using GWAS summary statistics can be written as Q B T S k = S k T W 1 M 2 , where W is a M × M diagonal matrix with W j j = 1 / M A F j 1 − M A F j , and M A F j is the minor allele frequency for the j t h variant (Lee et al., 2013). Under the null hypothesis, Q B T follows a scaled Chi-square distribution with one degree of freedom. In SKAT, the test statistic based on GWAS summary statistic is defined as Q S K A T S k = S k T WW T S k , where W is a M × M diagonal matrix of weights with W j j = B e t a M A F j ; a 1 , a 2 . The beta distribution density function has pre-specified parameters a 1 , a 2 and M A F j . 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 p k be the p-value of the k t h cohort study for k = 1 , … , K 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.

Suppose we have K cohorts with sample sizes n 1 , n 2 , … , n K , respectively. For the k t h cohort, the phenotype is modeled by a linear model y k = G k β k + ε k , where G k are matrices of genotypes with columns standardized to mean zero and variance 1, with dimension n k × M . y k is a n k × 1 vector of standardized phenotypes with mean zero and variance 1. β k is the vector of genotypes effect sizes in the specific gene and ε k is the vector of residuals representing environmental effects and non-additive genetic effects for the k t h 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 K cohorts, we suppose that the matrix of effect size β 1 , β 2 , … , β K has mean zero and covariance matrix V a r β 1 , β 2 , … , β K = 1 M   h 1 2 I ρ g 12 I … ρ g 1 K I ρ g 12 I … ρ g 1 K I h 2 2 I … ρ g 2 K I … … … ρ g 2 K I … h K 2 I , where ρ g i j is the genetic covariance on two phenotypes i and j on shared individuals, h k 2 is the heritability explained by the variants in a region for the k t h phenotype. The residuals ε 1 , ε 2 , … , ε K T has mean zero and covariance matrix V a r ε 1 , ε 2 , … , ε K T = 1 − h 1 2 I ρ e 12 I … ρ e 1 K I ρ e 12 I … ρ e 1 K I 1 − h 2 2 I … ρ e 2 K I … … … ρ e 2 K I … 1 − h K 2 I , where ρ e i j is the covariance of non-genetic effects on the i t h and j t h phenotypes among shared individuals. Then, the phenotypic correlation for cohort i and cohort j among the N s overlapping samples is ρ i j = ρ g i j + ρ e i j (Lemma). For these N s 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 β 1 , β 2 , … , β K and non-genetic component ε 1 , ε 2 , … , ε K T 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 y k = G k β k + ε k , k ∈ 1 , 2 , … , K .

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.

To evaluate the Type I error rates of Meta-TOW-S, we first simulate 1 , 000 regions, each encompassing approximately 100 genetic variants, reflecting the LD structure observed in individuals of European ancestry. We then replicate 10 5 times to generate the phenotypes under the null hypothesis of no genetic contribution to any of the three traits, that is β = 0 . Then we simulate 10 8 datasets to estimate the Type I error rates at nominal significance levels α = 2.5 * 10 − 6 , 10 − 5 , 10 − 4 . 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 1000 : 1000 : 1000 . In the second scenario, the sample sizes are unequal, with a ratio of 2000 : 1000 : 500 . For simplicity, the number of overlapping individuals between any two phenotypes is kept constant at 500 . The phenotypic correlation between any two phenotypes among the overlapping individuals is set to be either ρ = 0.5 or ρ = 1 . A correlation of 1 indicates simulation of mimicking one phenotype but form different cohorts

For Type I error evaluations, we use the significance levels α = 2.5 * 10 − 6 , 10 − 5 , 10 − 4 . 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.

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 α = 2.5 × 10 − 6 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 α = 2.5 × 10 − 6 . For the overlapping individuals, the phenotypic correlation between any pair of phenotypes across four scenarios: 0.1 , 0.2 , 0.5 , and 1 . 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 K = 3 , K = 5 and K = 10 . 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 h k 2 = 0.01 or h k 2 = 0.3 for k = 1 , … , K . The simulation results are depicted in Figures 1–4 for a heritability of 0.01 , while results for a heritability of 0.3 are presented in the Supplementary Material.

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 200 and genetic heritability h 2 = 0.01 . Specifically, we create datasets with 3 , 5 , and 10 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 1 , 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.

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 46 , 351 to 105 , 318 , 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 α = 2.5 × 10 − 6 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 w j = 1 / M A F j 1 − M A F j for Meta-Burden and w j = B e t a M A F j ; a 1 , a 2 for Meta-SKAT, where the beta distribution density function has pre-specified parameters a 1 , a 2 , and M A F j of the j t h 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.

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