Subset scanning for multi-trait analysis using GWAS summary statistics

Abstract Motivation Multi-trait analysis has been shown to have greater statistical power than single-trait analysis. Most of the existing multi-trait analysis methods only work with a limited number of traits and usually prioritize high statistical power over identifying relevant traits, which heavily rely on domain knowledge. Results To handle diseases and traits with obscure etiology, we developed TraitScan, a powerful and fast algorithm that identifies potential pleiotropic traits from a moderate or large number of traits (e.g. dozens to thousands) and tests the association between one genetic variant and the selected traits. TraitScan can handle either individual-level or summary-level GWAS data. We evaluated TraitScan using extensive simulations and found that it outperformed existing methods in terms of both testing power and trait selection when sparsity was low or modest. We then applied it to search for traits associated with Ewing Sarcoma, a rare bone tumor with peak onset in adolescence, among 754 traits in UK Biobank. Our analysis revealed a few promising traits worthy of further investigation, highlighting the use of TraitScan for more effective multi-trait analysis as biobanks emerge. We also extended TraitScan to search and test association with a polygenic risk score and genetically imputed gene expression. Availability and implementation Our algorithm is implemented in an R package “TraitScan” available at https://github.com/RuiCao34/TraitScan.


Introduction
Genome-wide association studies (GWAS) have successfully improved the understanding of the genetic basis of many traits.The emergence of deeply phenotyped GWAS databases such as UK Biobank (Bycroft et al. 2018), eMERGE (Gottesman et al. 2013), and Vanderbilt BioVU (Roden et al. 2008), has facilitated studying associations between single nucleotide polymorphisms (SNPs) and a large number of traits.Phenome-wide association studies (PheWAS) have utilized this rich source of SNP-trait relationships to explore disease risks (Denny et al. 2010) and drug development (Diogo et al. 2018).By evaluating each trait individually, PheWAS is computationally fast to implement.However, it is well documented that joint multivariate analyses can be more powerful than univariate analyses such as PheWAS (Robinson et al. 2018).Many efforts have been devoted to multi-trait analyses that evaluate the relationships between a SNP and a set of traits simultaneously (O'Reilly et al. 2012, Kim et al. 2015, Zhu et al. 2015, Li and Zhu 2017).Nonetheless, existing multi-trait analyses rely on domain knowledge to select a small number of related traits, and most of them only focus on obtaining high statistical power in hypothesis testing.
We are motivated to understand which risk factors contribute to childhood cancers, which are universally rare and have obscure disease etiology.For example, the etiology of Ewing sarcoma (EWS) remains unclear (Lahat et al. 2008), and conventional epidemiological methods to understand the disease are limited due to the extremely rare incidence rate of EWS (Spector et al. 2021).Recently, several genetic risk factors were identified in a GWAS of EWS (Machiela et al. 2018), potentiating the use of PheWAS to explore the risk factors agnostically.However, PheWAS can suffer from limited statistical power when scanning over a large number of traits, and simply applying the existing multi-trait methods to a large number of traits does not always yield meaningful results as the polygenic nature of some traits would eventually drive the statistical significance.To address the aforementioned limitations, we propose a novel multi-trait analysis method, TraitScan, that values both trait selection performance and high statistical power.Our method "TraitScan" is based on a fast subset scan framework (Neill 2012) with a linear scan time of the number of traits and thus can handle highdimensional trait selection.Our method contains three test statistics: higher criticism (HC), truncated chi-squared (TC), and a combined test of HC and TC.We note a similar method ASSET (Bhattacharjee et al. 2012) with the same objective.However, it requires an exhaustive search of all possible subsets, which results in an exponential scan time and thus is not computationally efficient when the number of traits exceeds a few dozen.
Our proposed TraitScan algorithm is able to utilize summary-level GWAS data in situations where individuallevel data are not available.Due to the logistical limitations and privacy concerns of sharing individual-level data, it has become a common practice to share GWAS summary statistics.Leveraging publicly available summary-level data, our method can filter relevant pleiotropic traits on any given SNP.Through simulations, we show that our method has high power and sensitivity in terms of trait selection under moderately sparse to sparse situations (i.e. the number of truly associated traits is smaller than or close to the square root of total traits).We evaluate traits associated with EWS through GWAS summary statistics on 754 traits filtered from the UK Biobank study (Sudlow et al. 2015).Besides single SNPs, we also show that our method can be extended to genetic scores, such as the predicted gene expression levels in transcriptomewide association studies (TWAS) (Feng et al. 2021) and polygenic risk scores (PRS) (Torkamani et al. 2018).We implement our method in an R package "TraitScan" that is publicly available at https://github.com/RuiCao34/TraitScan.The package provides an option to use the pre-calculated null distributions of the test statistics, which can handle screening 754 traits in 58 s.

Models
We start our formulation by considering continuous traits as outcomes of interest.Assume that the GWAS data consist of p continuous traits and a minor allele dose for a genetic variant of interest (i.e.SNP) collected from n individuals.For j ¼ 1; . . .; p, let y j ¼ ðy 1j ; . . .; y nj Þ T be a vector of values of the jth trait for the n individuals, x ¼ ðx 1 ; . . .; x n Þ T a vector of the minor allele doses of the SNP of interest, and j ¼ ð 1j ; . . .; nj Þ T a vector of error terms for the jth trait.We assume that each trait can be modeled as a linear function of the genetic variant, and without loss of generality, y j is centered to have mean 0: (1) Our method performs a scan for traits under the null hypothesis that none of the traits is associated with the SNP: and the alternative hypothesis that at least one trait is associated with the SNP: (3) where j 2 f1; . . .
We can stack the models together as: To model the correlation structure among the continuous traits, we assume that x is fixed and the rows of represent n Â p i.i.d.observations from an MVNð0; XÞ distribution, where 0 is a vector of zeroes of length n Â p and X ¼ I nÂn R, where R is the potentially unknown p Â p covariance matrix of the traits, i.e. ð i1 ; . . .; ip Þ $ MVNð0; RÞ for an arbitrary individual i.When the individual-level data are available, the matrix can be estimated from the residuals of fitting separate linear regression models: R ¼ 1 nÀ1 ðY À x bÞ T ðY À x bÞ, where b is the ordinary least square estimate.
In many cases, only summary statistics are available, i.e. instead of having individual-level data, we have access to bj , seð bj Þ from Equation ( 1) and the z-score from a Wald test z j ¼ bj =seð bj Þ.We approximate the ðj 1 ; j 2 Þ th entry of R, rj 1 j 2 by using the null SNPs (i.e. the SNPs with no association with any traits) and ignoring the estimation error of b (Kim et al. 2015, Liu andLin 2018): In addition, when the summary statistics come from overlapping but not identical samples, the correlation between z-scores is still proportional to the trait correlation r j1j2 (Li et al. 2021): where n j 1 ; n j 2 , and n j 1 j 2 are the sample sizes of trait j 1 , trait j 2 , and their overlapping samples respectively.Since the ðj 1 ; j 2 Þ th z-score correlation is a constant across all null SNPs, it can be estimated empirically (Zhu et al. 2015): where z k j 1 is the z-score for null SNP k and trait j 1 , and z j 1 is the mean of vector ðz 1 j1 ; ::z k j1 ; . ..Þ across all null SNPs.For binary traits, the statistical model and z-score correlation estimation can also be generalized.In GWAS, the logistic regression models are usually fitted as: When the effect size b j1 is small, which most often happens in GWAS, the logistic regression model can be approximated by a linear regression model based on the first-order Taylor expansion on b j1 : where a j0 , a j1 , and n can be regarded as linear regression coefficients.The covariance for binary traits can be similarly derived using summary statistics of the null SNPs.

Scan statistics
To efficiently detect the most anomalous subset of traits among a large number of putative traits, we use the linear time subset scan (LTSS) framework (Neill 2012), which only requires searching a linear number of traits (p) rather than an exponential number of traits (2 p ) across all possible trait combinations.The LTSS framework defines a score function FðSÞ that measures the anomalousness of the trait subset S with a corresponding priority function G.Under the strong LTSS property, to maximize FðSÞ, we could simply sort traits by G and compare FðSÞ by different numbers of traits (from 1 to p, details in Supplementary Section S1.1).We accordingly construct two score functions, higher criticism (HC) and truncated chi-squared (TC).We show by simulations that the HC or TC method had better performance than PheWAS and other competing methods under different scenarios.We also combine the two tests by taking the minimum P-value of the two tests, which allows us to attain results comparable to the better-performed HC or TC method in terms of statistical power and trait selectivity.The overview of TraitScan algorithm is present in Fig. 1.

Decorrelation
As in the LTSS framework, a trait-level statistic is required to quantify the amount of association between a trait and a genetic variant.We used the P-value p j of the Wald test z j from the separate regression models [Equation ( 1)], for j ¼ 1; . . .; p.As the traits are correlated and sampled from the same or overlapping individuals in our framework, the p j 's are correlated.Herein, we perform the ZCA-cor whitening method (Kessy et al. 2018) on z j , ensuring that the whitened z Ã j remains maximally correlated with z j .Then we obtain p Ã j corresponding to z Ã j .aÞ and jjSjj is the cardinality of the trait subset S. The score function is the standardized difference between the observed count of P-values lower than a P-value threshold a and the expected count.According to Theorem 1 (Supplementary Section S1.1), it can be proved that the HC statistic satisfies the strong LTSS property with priority function G HC;a ðp Ã j Þ ¼ Iðp Ã j < aÞ.As we do not know the optimal a, we define grid-based HC test statistic H HC :

Higher criticism statistic
over a grid of a and its corresponding subset An ideal a grid should ensure that all possible subsets are in the search space, i.e. there should be no more than one P-value between two arbitrary adjacent a's.In practice, we recommend using the a grid as a geometric sequence from the Bonferroni significant P-value threshold to overall Type I error with a sequence length of 200, i.e. a 1 ; . . .; a 200 ¼ 0:05=p; . . .; 0:05 |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} length¼200 .The lower bound of a ensures that our test would always be more powerful than PheWAS in special scenarios where all traits are uncorrelated.In the meantime, we use an upper bound of 0.05 to decrease the search space.

Truncated chi-squared statistic
The HC test may not have ideal performance under nonsparse scenarios (Barnett et al. 2017) and does not take into account the strength of association.To make our method more robust, we propose an additional statistic, which is similar to the truncated z-score method (Bu et al. 2020) and also meets the strong LTSS property.First, we define c as the jz Ã j threshold, which is closely related to a in HC statistics: c ¼ U À1 ð1 À a=2Þ and U as the cumulative distribution function of a standard normal distribution.The score function FðSÞ is defined as F TC;c ðSÞ ¼ ÀlogðP McðSÞ jH 0 Þ, i.e. the negative log Pvalue of the subset score function M c ðSÞ ¼ The TC statistic also meets the condition of Theorem 1 (Supplementary Section S1.1) and thus satisfies the strong LTSS property.
Similarly, we test a grid of c one-to-one mapping to the grid of a defined previously: c grid is chosen in correspondence with the a values.

Assessing significance
We have now obtained two subsets S HC and S TC that maximize HC or TC statistics.As described above, both subsets always include at least one trait.To determine whether the selected subset is sufficiently anomalous, we calculate the corresponding P-values p HC and p TC by comparing the two statistics H HC and H TC with their distributions under the null hypothesis.A trait scanning method p HC can be calculated analytically (Supplementary Section S1.2).As for p TC , we use Monte Carlo (MC) simulations.We start by simulating p z-scores under the null, i.e. standard normal distribution for B iterations.For the bth iteration, the H TC;b can be calculated, and empirical P-values p TC can be estimated from the simulated null distribution.
To combine the HC and TC tests, we compare the P-values p HC and p TC and get the grid-based statistics: and the traits selected by the combined test are also determined by the test with a smaller P-value: where M ¼ HC or TC.The empirical null distribution of H combined can be similarly simulated by MC, and the P-value from the combined test p combined is calculated by comparing the test statistic H combined and its distribution under the null.
If the null hypothesis is rejected, we can conclude that S combined is the most anomalous trait subset and the SNP is associated with at least one of the traits in S combined .Note that this MC simulation step only depends on the number of traits p and the choice of a grid.For SNPs sharing the same number of traits p, we do not need to recompute the test statistic distribution under the null.

Extension to genetic scores
Genetic scores integrate information from multiple SNPs.Linking traits with genetic scores could bring in more statistical power and provide a meaningful interpretation of the results.Genetic scores, which are usually the linear combinations of allele counts of multiple SNPs, have been extensively developed and distributed.Polygenic risk scores that predict the risk of clinical and epidemiological traits (Lewis and Vassos 2020) or imputation models for gene expression levels in TWAS (Xu et al. 2023) are two types of commonly used genetic scores.We will show how TraitScan can be easily utilized on genetic scores using summary-level GWAS data and an external genetic reference panel.

Continuous traits
Let X gs denote the genetic score from q SNPs: X gs ¼ P q l¼1 c l X l , where X l ¼ ðx 1l ; . . .; x nl Þ is the genotype vector for n independent individuals at the l th SNP, and c l is the SNP weight vector.In GWAS models, the j th trait y j is marginally regressed on each SNP X l , and regression coefficients bjl and seð bjl Þ are estimated from the linear regression model For the genetic score, we are interested in the regression model (20) where r2 gs is the residual variance estimate and n j is the sample size for jth trait.The items X T gs y j and y T j y j can be derived from GWAS summary data (Pattee and Pan 2020): where ŝ2 l is the variance of SNP l.Both ŝ2 l and the genotype matrix X T gs X gs can be estimated from a reference panel comprising genotypic data of individuals from a general population (Auton et al. 2015).In practice, for Equation ( 23), we can calculate y T j y j across multiple SNPs and take the median as the estimate.
After z statistics fz j;gs g are computed, we could follow the same decorrelation and trait scanning steps as above since the genetic score can also be treated as a SNP, and the covariances between fz j;gs g are identical under the null hypothesis.

Binary traits
We have the logistic regression for binary traits: for ith individual, jth trait, and lth SNP, b 0;jl and b jl are the regression coefficients.Following Pattee and Pan (2020), we could approximate Pðy ij ¼ 1jx il Þ as a continuous outcome under a linear regression model and denote b 0;jl and b jl as the coefficients.The following equations hold: bjl ¼ e À b0;jl ð1 þ e À b0;jl Þ 2 bjl ; (25) where e À b0;jl ¼

Pðyij¼0Þ
Pðyij¼1Þ is the ratio of control and case sizes.The logistic regression coefficients can be thus converted to linear regression coefficients and handled by the steps mentioned above.
We used our method on UK Biobank GWAS data to find out potential traits linked to EWS.EWS is a type of rare childhood cancer in bone or soft tissue (Li and Chen 2022).Previous studies (Postel-Vinay et al. 2012, Machiela et al. 2018) suggested that six SNPs rs113663169, rs7742053, rs10822056, rs2412476, rs6047482, and rs6106336 were significantly associated with EWS in individuals of European ancestry.We analyzed these six SNPs using the GWAS summary statistics of the UK Biobank data.
UK Biobank is a large-scale database encompassing a broad range of phenotypes, where individuals' genetic data are linked to electronic health records and survey measures (Sudlow et al. 2015).More method and analysis details can be found on the Pan-UK Biobank website (https://pan.ukbb.broadinstitute.org/).The UK biobank GWAS summary-level data analyzed in this study was downloaded after the 3 August 2023 update.
Since EWS occurs in both males and females of European ancestry, we focused on the GWAS summary statistics for individuals of European ancestry and with a sufficient number of participants of both sexes.We applied the following criteria which left us with 754 traits to perform the TraitScan algorithm: • Continuous traits with a sample size of at least 5000 or binary traits with a sample size of at least 5000 cases and 5000 controls.• Traits with genetic heritability P-value < 0:05 by stratified LD score regression (S-LDSC) (Finucane et al. 2015).
• Traits with the sample size of each sex larger than 20.
• Traits belonging to these categories were included: healthrelated outcomes, online follow-up, biological samples, X-ray absorptiometry (DXA), cognitive function, verbal interview, touchscreen questions (except traits related to eyes), and baseline characteristics.
As we intended to evaluate six SNPs, the significance level for TraitScan tests was set at 0:05=6 ¼ 0:0083 after the Bonferroni correction.We compared TraitScan with PheWAS with Bonferroni adjustment (significance level 0:05=ð754 Â 6Þ ¼ 1:11 Â 10 À5 ).Given the large number of putative traits, no other competing multi-trait analysis method was applied due to the prohibitive computational time.As a result, SNPs rs113663169, rs10822056, rs2412476, rs6047482, and rs6106336 were shown to have significant associations with at least one trait out of 754 examined traits in UKBiobank by TraitScan (TraitScan combined test P-values 1 Â 10 À4 ).For SNP rs7742053, TraitScan combined test P-value was 0.872 and thus did not reach statistical significance.Table 1 summarizes the results of TraitScan combined test for the most significant trait in each category identified for the five SNPs.It showed traits that were highly significant in PheWAS were also captured by TraitScan.In fact, TraitScan identified a total of 28 UK Biobank traits related to the five EWS-linked SNPs, while 8 of the trait-SNP associations did not reach statistical significance in PheWAS.A full list of selected traits is shown in Supplementary Materials (Supplementary Table S1).
To demonstrate the TraitScan application on genetic scores, we further carried out our method on the transcriptomic scores of gene KIZ and gene RREB1, i.e. the SNP-imputed gene expression of the two genes, with the 754 UK Biobank traits.KIZ and RREB1 had strong evidence to be linked to three top genome-wide significant SNPs in EWS GWAS (KIZ was linked to SNPs rs6047482 and rs6106336, and RREB1 to rs7742053) (Machiela et al. 2018).The weights in the transcriptomic scores of human blood were obtained from Xu et al. (2023), and the internal r 2 of the scores were 0.206 and 0.031 for KIZ and RREB1, respectively.The significance levels for both genes were set at 0:05=2 ¼ 0:025 after the Bonferroni correction.In addition, we also tested the relationship between the number of risk alleles identified in Machiela et al. (2018) with EWS.It was shown that EWS cases had on average 1.08 more risk alleles than controls (P-value ¼ 2:44 Â 10 À63 ).
After applying TraitScan on the genetic scores, neither of the two genes KIZ (P-value ¼ 0:0542) and RREB1 (P-value ¼ 1) reached statistical significance in TraitScan, although the trait insulin-like growth factor 1 (IGF-1) picked up by KIZ was marginally significant.On the other hand, imputed gene expression of RREB1, of which rs7742053 was an expression quantitative trait loci (eQTL), had no evidence of linking to any of the 754 traits.The genetic score of six EWS SNPs, however, was significantly associated with three traits: hair color (blonde), skin color, and facial aging (TraitScan Pvalue < 1 Â 10 À4 ), while the PheWAS identified two additional traits on the EWS score: hair color (dark brown) and ease of skin tanning.
To investigate the causal relationship between EWS and the selected traits, we further performed bidirectional Mendelian A trait scanning method randomization (MR) analysis on the traits selected by TraitScan using the TwoSampleMR package (Hemani et al. 2017(Hemani et al. , 2018)).The instrumental variables were selected from either UK Biobank or EWS GWAS data and were clumped by r 2 < 0:001 and P-value < 5 Â 10 À5 .For the selected traits with more than one instrumental variable, inverse-variance weighted (IVW), Egger regression, weighted median, simple mode, and weighted mode methods were applied, while the Wald ratio method was applied for the selected traits with one instrumental variable.Full results are reported in Supplementary Table S2.Lower alkaline phosphatase levels, the trait selected to be associated with SNP rs10822056 in TraitScan but not PheWAS, was shown to be causal for EWS (MR weighted mode coefficient ¼ À1:22 with P-value ¼ 2:07 Â 10 À4 ).For the causal direction where EWS was the exposure, no MR test reached statistical significance after accounting for multiple testing, suggesting no causal effect from EWS to the selected UK Biobank traits.

Simulation
Throughout the simulations, we used five metrics to assess the performance of each method, i.e. power, size, recall, precision, and Jaccard similarity.Let p be the total number of traits, S Ã be the subset of traits as chosen by a particular method and let S 0 be the true subset of pleiotropic traits.The size was defined as jjS Ã jj, which is the number of selected traits by a multi-trait method.Then, we defined precision to be jjS Ã \S0jj jjS Ã jj , the proportion of traits identified by the method that was truly associated with the SNP.We defined recall to be jjS Ã \S0jj jjS 0 jj , the proportion of the pleiotropic traits identified by the method.We defined Jaccard similarity, a combination of precision and recall, to be jjS Ã \S0jj jjS Ã [S0jj .Finally, power was assessed by comparing the observed statistic to the simulated distribution of null statistics.We reported the average value of all the metrics across the simulation iterations.
We compared the performance of variable selection and testing for TraitScan using summary statistics with some existing methods: PheWAS (with Bonferroni adjustment), CPASSOC (S hom , S het ) (Zhu et al. 2015), MTaSPUs (Kim et al. 2015), and generalized higher criticism (GHC) (Barnett et al. 2017).Among these methods, MTaSPUs and S hom cannot select traits, and thus we only evaluated their performance in terms of statistical power.S het includes a parameter grid as the thresholds for z-scores and can naturally select out the traits with absolute z-scores smaller than each threshold.As suggested by their package, the parameter grid was set as the observed trait P-values.The GHC was originally proposed for SNP set association with a single trait, and here we used it for a single SNP association with multiple traits.We did not compare our method with ASSET because it is not computationally efficient in our simulation settings (754 traits for the real data variance-covariance scenario and 50 traits for the rest of the scenarios).
We conducted simulations under multiple scenarios to show TraitScan has higher power and Jaccard similarity under moderately sparse (jjS 0 jj=p < 0:4) and sparse situations (jjS 0 jj=p < 0:05).We focus the discussion on the scenarios with varying numbers of truly associated traits (scenario 1) or with real data variance-covariance matrix (scenario 2) and briefly discuss the other four scenarios and their results.
We assessed the method performance by varying the number of truly associated traits jjS 0 jj in scenario 1 (Fig. 2).In terms of statistical power, we observed that TraitScan test statistics (HC, TC, HCþTC) had the highest power under moderately sparse and sparse situations (jjS 0 jj < 22, or jjS 0 jj=p < 0:44).The HC test was more powerful than the TC test under extremely sparse situations (jjS 0 jj ¼ 1).When jjS 0 jj ¼ 50, the GHC, MTaSPUs, PheWAS, and S hom had the highest power, and TraitScan and S het were less powerful.In terms of variable selection performance, we found that TraitScan test statistics (HC, TC, HCþTC) also had the highest Jaccard similarity under moderately sparse and sparse situations (jjS 0 jj < 22, or jjS 0 jj=p < 0:44), while S het had better Jaccard similarity as increasing proportion of truly associated traits.However, S het tended to over-select traits and thus had the lowest precision under most situations.
We also tested method performance using the covariance matrix and effect sizes estimated from the 754 UK biobank traits (Fig. 3), where 74 traits with marginal P-value 0:05 were set as a true subset of pleiotropic traits.In this scenario, we compared the performance of PheWAS and TraitScan under trait-SNP associations of different strengths.S het , GHC, and MTaSPUs were not applied due to the computational burden.In addition, S hom was excluded because the simulation setting evidently did not adhere to the homogeneous effect assumption, an assumption that is rarely met in real data analysis.By comparing the Type I errors (b ¼ 0), we showed that TraitScan tests were well calibrated.Similar to scenario 1 mentioned above, TraitScan had higher power and Jaccard similarity than PheWAS under all effect sizes.We also noticed that the selected trait size did not grow with effect size.When the effect sizes were small, the P-values of truly associated traits and null traits were close, and TraitScan tended to select a large set of traits as the most anomalous trait subset.
We showed that TraitScan can handle mixed types of traits (i.e.continuous and binary) simultaneously (Supplementary Fig. S1) and the performance followed the same pattern as continuous traits.TraitScan still demonstrated the highest power and Jaccard Similarity under effects varying in both directions and magnitudes (Supplementary Table S4), under block-diagonal correlated traits (Supplementary Fig. S2), or different correlation magnitudes or structures (Supplementary Table S5).Moreover, we found that TraitScan had higher power and Jaccard similarity when handling more highly correlated traits.The detailed parameter settings of the simulations can be found in Supplementary Section S3.1.

Discussion
We proposed a new method called TraitScan for post-GWAS trait subset scanning and testing.While most of the existing multi-trait methods rely on domain knowledge, our method allows agnostic search among a large number of traits and is able to identify a set of traits with the most anomalousness.TraitScan utilizes the fast subset scan framework (Neill 2012), resulting in a linear scan time over the number of traits.Taking correlation among traits into consideration, TraitScan has demonstrated higher power and trait selectivity than PheWAS when sparsity was low or modest.The method is compatible with both individual-level and summary-level GWAS data, although we focus more on summary-level GWAS data herein to allow an easy application to existing deeply phenotyped GWAS summary statistics databases.
Given a fixed number of a 0 s, TraitScan has an OðpÞ time complexity.If the number of a 0 s is proportional to p, the time complexity is Oðp 2 Þ.While for the other multi-trait analysis methods, the S het test in CPASSOC also has an OðpÞ time complexity given a fixed P-value threshold, since it ranks the P-values and directly applies the threshold onto the raw P-values.When the thresholds are set as the input P-values, which is recommended in the CPASSOC pipeline, its time complexity is also Oðp 2 Þ.The ASSET has an Oð2 p Þ time complexity, meaning the computational time will be doubled once adding one more trait.Table 2 lists the computational time for analyzing 10 traits 5000 times using different methods, including TraitScan with 10 000 iterations in the MC simulation (TraitScan-MC), TraitScan: HC test based on the analytical null distribution (TraitScan-analytic), and TraitScan test using a precalculated null distribution estimated by MC with a given p (TraitScan-precalculated). All the programs  were ran on Intel Haswell E5-2680v3 processors with 16-GB memory, and the trait correlations and b parameters were the same as in Scenario 1.
In the examination of traits associated with EWS, TraitScan identified eight additional trait-SNP associations which did not reach the PheWAS significance level.One of these traits, alkaline phosphatase, measured by blood assays, also showed significance in MR analysis, suggesting it was causally related to EWS.Evidence has shown the presence of abundant alkaline phosphatase activity in EWS tumor cells (Sharada et al. 2006); however, the direction of association between alkaline phosphatase and EWS was previously unknown.Another trait, IGF-1, was selected by TraitScan for rs6047482 and rs6106336 on chromosome 20.IGF1-receptor is known to be upregulated in EWS, and anti-IGF1 is an experimental therapy (Gonzalez et al. 2020).Besides, the SNP rs7742053, the only SNP that failed to reach genome-wide significance in both TraitScan and PheWAS, has recently been reported to have a specific role in the increased binding of GGAA microsatellite alleles with the chromosomal translocation encoding chimeric transcription factors (Lee et al. 2023).
Examining the genetic score of KIZ and RREB1 allowed us to investigate whether any trait was associated with EWS on the gene level.If a gene is associated with the same set of traits, then likely multiple SNPs in the gene will be associated with the traits, leading to higher power than the single SNP test.However, if the weight in the genetic score is not informative, such as imputing gene expression in a non-relevant tissue, or if multiple SNPs in the gene suggest a different association with the trait, we would have diminished power.We did not observe any traits significantly associated with the genetic scores (i.e.imputed gene expression) of KIZ and RREB1, potentially because blood may not be the most relevant tissue for EWS.We note that like other methods, our results relied on the quality of GWAS data.When handling real GWAS data, we applied a couple of filtering steps to exclude traits that are not heritable.However, after the filtering steps, there were still a few traits that lacked reasonable explanations of their genetic heritability such as the inpatient record format, or the potential mechanisms to be associated with EWS, such as fruit intake within the past 24 h.We suspect it was due to the inadequate adjustment for confounding in the original GWAS analysis (Holmes et al. 2019).Without accessing the individual-level data, it is difficult to examine or correct the summary-level GWAS data, although there is some recent work performing quality control on GWAS errors using summary statistics and a reference panel (Chen et al. 2021, Darrous et al. 2021).
To use TraitScan in real data analysis, the following additional steps could help avoid potential power loss and increase the interpretability of the results.First of all, we suggest removing highly correlated traits from the pool of putative traits by examining the empirical trait correlation matrix.The strong LTSS property of TraitScan requires traits to be independent of each other.As shown in the simulation, our method had relatively low statistical power when the genetic variant had the effects and correlations of the same direction on most of the traits.We find that the decorrelation step on z-scores shifted the means of z-scores of the truly associated traits toward zero, resulting in a power loss.Therefore, removing such traits could potentially improve the statistical power.Future work may be focused on developing subset algorithms balancing the computational time and scan sensitivity.Secondly, we recommend checking the correlation between the decorrelated traits and raw traits.Due to the trait decorrelation in TraitScan, trait selection and testing are performed on the decorrelated z-scores, which are essentially linear combinations of raw z-scores.Although the ZCA-cor decorrelation method maximizes the average correlation between each dimension of the decorrelated and original data, the decorrelated traits might be considered to differ from the original traits.Therefore, this step could improve the interpretability of the findings.In our real data analysis, 99% of the 754 UK Biobank traits had an empirical correlation with the original trait greater than 0.7.
The understanding of rare diseases such as childhood cancer has long been limited.TraitScan is able to provide a list of possible traits associated with EWS through the diseaselinked genetic variants.As association does not imply causality, further biological experiments or additional data analysis approaches such as MR are required to study whether a trait and target disease are causally linked and whether the trait is a risk factor or a consequence of the disease.
When individual-level data are available, the regression coefficients bj;gs and seð bj;gs Þ with z statistic from the Wald test z j;gs ¼ bj;gs =seð bj;gs Þ can be directly calculated.When only summary-level data are available, we have 17) and test the null hypothesis H 0;gs : b 1;gs ; . . .; b p;gs ¼ 0: (18)

Table 1 .
Use TraitScan to search among 754 traits in UK Biobank for EWS-linked SNPs.