Joint Analysis of Multiple Traits Using "Optimal" Maximum Heritability Test

The joint analysis of multiple traits has recently become popular since it can increase statistical power to detect genetic variants and there is increasing evidence showing that pleiotropy is a widespread phenomenon in complex diseases. Currently, most of existing methods use all of the traits for testing the association between multiple traits and a single variant. However, those methods for association studies may lose power in the presence of a large number of noise traits. In this paper, we propose an “optimal” maximum heritability test (MHT-O) to test the association between multiple traits and a single variant. MHT-O includes a procedure of deleting traits that have weak or no association with the variant. Using extensive simulation studies, we compare the performance of MHT-O with MHT, Trait-based Association Test uses Extended Simes procedure (TATES), SUM_SCORE and MANOVA. Our results show that, in all of the simulation scenarios, MHT-O is either the most powerful test or comparable to the most powerful test among the five tests we compared.


Introduction
Increasing evidence shows that pleiotropy, the effect of one variant on multiple traits, is a widespread phenomenon in complex diseases [1]. Furthermore, in genetic association studies of complex diseases, multiple related traits are usually measured. For example, hyperuricemia is usually present in patients with gout [2]; coronary heart disease is predicted by cytokine interleukin-6, C-reactive protein, interleukin-1, tumor necrosis factor-α and fibrinogen [3,4]; and neuropsychiatric disorders depend on a range of overlapping clinical characteristics [5]. Although most published genome-wide association studies (GWASs) analyze each of the related traits separately, joint analysis of multiple traits may increase statistical power to detect genetic variants [6][7][8][9]. Thus, joint analysis of multiple traits has recently become popular.
Several statistical methods have been developed for joint analysis of multiple traits. These methods can be roughly divided into three groups: combining the univariate analysis results, regression methods, and dimension reduction methods. For combining univariate analysis results, one first conducts the univariate test by performing an association test for each trait individually and then combines the univariate test statistics or combines the p-values of the univariate tests [2,[10][11][12]. Regression methods include mixed effect models [9,13,14], generalized estimating equation (GEE) methods [15,16], and reverse regression methods [5,17]. Mixed effect models can account for relatedness, population structure, and polygenic background effect, but it is computationally challenging. The GEE methods, based on a marginal regression model, allow the variant having different effect sizes and effect directions on different traits. These methods can also accommodate covariates and different types of traits. Reverse regression methods take genotypes as the response variable and multiple traits as independent predictors, therefore, reverse regression models do not need to know the complex distributions of traits and can be applied to a large number of mixed types of traits. Dimension reduction methods include canonical correlation analysis (CCA) [18], principal components of traits (PCT) [19], and principal components of heritability (PCH) [20][21][22][23]. CCA is to seek a linear combination of multiple variants and a linear combination of multiple traits such that the correlation between the two linear combinations reaches its maximum. The PCT methods are usually based on the first PC or first few PCs of the traits [22,24]. However, as Aschard et al. [2014] showed that testing only the first few PCs often has low power, whereas combining signals across all PCs can have greater power. Nevertheless, it is not clear how many PCs are needed, and how robust these methods are when there exists noise traits. PCH is to find a linear combination of multiple traits such that this linear combination has the maximum heritability.
In this article, we first propose a maximum heritability test (MHT). Based on MHT, we develop an "optimal" maximum heritability test (MHT-O) to test the association between multiple traits and a single variant. In each step of MHT-O, we delete one trait that has the weakest association with the variant. Then, we find the optimal number of traits and use MHT to test the association between the optimal number of traits and the variant. Using extensive simulation studies, we compare the performance of MHT-O with MHT, Trait-based Association Test uses Extended Simes procedure (TATES) [11], SUM_SCORE and MANOVA [8]. Our results show that, in all of the simulation scenarios, MHT-O is either the most powerful test or comparable to the most powerful test among the five tests we compared.

Method
We consider a sample with n unrelated individuals. Each individual has K (potentially correlated) traits and has been genotyped at one variant. Let Y = (Y 1 ,. . .,Y K ) T denote the random vector of K traits and X denote the random variable of the genotype score at a variant. Let y i = (y i1 ,. . .,y iK ) T denote the values of K traits and x i denote the genotype score of the i th individual, where x i is the number of minor alleles that the i th individual has at the variant. We can consider that y 1 ,. . .,y n is a random sample from Y and x 1 ,. . .,x n is a random sample from X. Now, let us consider linear models We partition the total phenotypic covariance of Y as V P = V G + V R [25]; V G = var[β 1 X,. . ., β K X] = var(X)ββ T is the genetic variance due to the genotype scores X, where β = (β 1 ,. . ., β K ) T ; V R = var[ε 1 ,. . ., ε K ] is the residual covariance after removing the genetic effect. var(X) can be estimated by 1 x i . β and V R can be estimated from the linear models β k is estimated by the least square estimator. Let r ik denote the estimates of residuals ε ik . Then, the (j, k) th element of V R is estimated by 1 n X n i¼1 r ij r ik .
Let us consider a linear combination of Y, w T Y ¼ heritability of w T Y can be written as The heritability of w T Y depends on w and we can find a linear combination of w T Y that has the largest heritability among all linear combinations of Y. We define the maximum heritability as the test statistic to test the association between these K traits and the variant. We denote this test as maximum heritability test (MHT). The MHT statistic can be written as where λ max (A) denotes the largest eigenvalue of matrix A. However, the test statistic T MHT may lose power in the presence of a large number of noise traits. Therefore, we propose an "optimal" maximum heritability test (MHT-O) to test the association between multiple traits and the variant. MHT-O includes a procedure of deleting traits that have weak or no association with the variant. It has the following steps: Step 1. Given traits Y = (Y 1 ,. . .,Y K ), initialize r = K and Y (r) = Y. Denote T MHT, r as T MHT based on Y (r) .
Step 2. Denote T Ài MHT ; r as T MHT based on Y (r) with the i th trait deleted for i = 1,. . .,r; denote I ¼ arg max i T Ài MHT ; r and T MHT ; rÀ1 ¼ T ÀI MHT ; r . Let Y (r−1) denote Y (r) with the I th trait deleted and update r = r − 1.
Denote p r as the p-value of T MHT, r . The test statistic of MHT-O is defined as We use a permutation test to evaluate the p-value of T MHT−O . Intuitively, two layers of permutations are needed to estimate p r and the overall p-value for the test statistic T MHT−O . Ge et al. [26] proposed that one layer of permutation can be used to estimate these p-values. We use the permutation procedure of Ge et al. to estimate p r and the overall p-value for the test statistic T MHT−O . In details, we randomly shuffle the genotypes in each permutation. Suppose we perform B times of permutations. Let T ðbÞ MHT; r denote the value of T MHT, r based on the b th permuted data, where b = 0 represents the original data. Then, we transfer T ðbÞ MHT; r to p ðbÞ r by

Comparisons of Methods
We compare our proposed method with MHT, TATES [11], MANOVA [8], and SUM_SCORE. TATES combines p-values obtained in a standard univariate GAWS to acquire one trait-based pvalue, while correcting for correlations between components. SUM_SCORE performs an association test for each trait individually to obtain the univariate score test statistic for each trait. Then, the test statistic of SUM_SOCRE is the summation of the univariate score test statistics. We use asymptotic distributions to evaluate the p-values of SUM_SCORE, TATES and MANOVA.

Simulation
To evaluate the type I error rates and powers of MHT and MHT-O, we generate genotypes according to minor allele frequency (MAF) and assume Hardy Weinberg equilibrium. Then, we generate K traits by the factor model [11,19] where y = (y 1 ,. . .,y K ) T ; x is the genotype score at the variant of interest; λ = (λ 1 ,. . .,λ K ) is the vector of effect sizes of the genetic variant on the K traits; f = (f 1 ,. . .,f R ) T * MVN(0, S), S = (1 − ρ) I + ρA, A is a matrix with elements of 1, I is the identity matrix, and ρ is the correlation between factors; γ is a K by R matrix; c is a constant number; and ε = (ε 1 ,. . ., ε K ) T is a vector of residuals, and ε 1 ,. . ., ε K are independent, and ε k * N(0, 1) for k = 1,. . ., K.
To evaluate type I error rates of MHT and MHT-O, we let β = 0. To evaluate powers, we let β > 0. In the simulation studies for evaluation of type I error rates and powers, we set MAF = 0.3 and ρ = 0.2.

Results
To evaluate the type I error rates of the two proposed tests (MHT and MHT-O), we consider 20 quantitative traits. We also consider different sample sizes, different significance levels, and different models. In each simulation scenario, the p-values of MHT and MHT-O are estimated by 1,000 permutations and the type I error rates of the two tests are evaluated using 10,000 replicated samples. For 10,000 replicated samples, the 95% confidence intervals (CIs) for estimated type I error rates of nominal levels 0.05 and 0.01 are (0.046, 0.054) and (0.008, 0.012), respectively (see Appendix for details). The estimated type I error rates of the two tests are summarized in Table 1. From this table, we can see that 58 out of 60 (greater than 95%) estimated type I error rates are within the 95% CIs and the two estimated type I error rates (0.05415 and 0.0126) not within the 95% CIs are very close to the bound of the corresponding 95% CI, which indicates that the two tests are all valid.    (Figs 2 and 3) are similar to that for 20 traits (Fig 1). We also give power comparisons of the five tests using a significance level of 5×10 −8 with 10 8 permutations and 500 replicates for 20 traits under model 1 (S1 Fig). S1 Fig  shows that the patterns of the power comparisons using significance level 5×10 −8 are similar to that using a significance level of 0.05 in Fig 1 (model 1). In summary, MHT-O is either the most powerful test or comparable to the most powerful test among all the tests we compared. Therefore, our MHT-O is a robust test to a variety of models.

Discussion
We propose MHT-O to perform joint analysis of multiple traits in association studies based on the following reasons: (1) multiple related traits are usually measured in genetic association studies of complex diseases; (2) there is increasing evidence showing that pleiotropy is a widespread phenomenon in complex diseases; and (3) the power of existing methods decreases in the presence of non-associated traits. The proposed MHT-O includes a procedure of deleting traits that have weak or no association with the variant. Therefore, it can be robust to the existence and the number of non-associated traits. By deleting one trait that has the weakest association with the variant in each step, MHT-O can maintain high power in the presence of a large number of non-associated traits. This feature is essentially important when there exist a large number of correlated traits but there are no guidelines to select relevant traits. Our results show that MHT-O has correct type I error rates and is either the most powerful test or comparable to the most powerful test among the five tests we compared. No other methods in the simulation studies show consistent good performance.
Due to the allelic heterogeneity and the extreme rarity of individual variants in rare variant association studies, the variant-by-variant methods for common variant association studies may not be optimal [27]. It has been shown by recent studies that complex diseases are caused by both common and rare variants [28][29][30][31][32][33][34]. Statistical methods including burden tests [27,[35][36][37][38], quadratic tests [39][40][41], and combined tests [42][43][44] have been developed for rare variant association studies with a single trait. Currently, there are limited researches on rare variant association studies for joint analysis of multiple traits [14,45]. MHT-O can be extended to rare variant association studies by extending Eq (1) to include multiple variants. MHT-O can also be extended to family-based studies by extending Eq (1) to mixed linear model. However, the performance of MHT-O in rare variant association studies and in family-based association studies needs further investigation.
The fact that population stratification can seriously confound association results has been long recognized in association studies based on unrelated individuals [46,47]. Several methods to control for population stratification have been developed for association studies based on unrelated individuals. These methods include principal component (PC) approach [48][49][50][51][52], genomic control (GC) approach [53][54][55], and mixed linear model (MLM) approach [29,56]. Like most association tests based on unrelated individuals, MHT-O subjects to bias due to population stratification. To make MHT-O robust to population stratification, we can use the PC approach. Let P i = (p i1 ,. . .,p iL ) T denote the first L PCs of the genotypes at a set of genomic markers for the i th individual. Let y Ã ik and x Ã i denote the residuals of the regressions y ik ¼ a 0k þ a T k P i þ ε ik and the residuals of the regression x i = α 0 + α T P i + ε i , respectively. Using y Ã ik and x Ã i to replace y ik and x i , we can make MHT-O robust to population stratification. However, the performance of using the PC approach to control for population stratification in MHT-O needs further investigations. where α is the significance level. Then, Pr(ξ = 1) = α and Pr(ξ = 0) = 1 − α because p follows a uniform distribution between 0 and 1 under the null hypothesis. Suppose there are R replicates. Let ξ i denote the value of ξ for the i th replicate, where i = 1,. . .,R Then, the estimated type I error rate is given by x We define ða À 1:96 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 À aÞ=R p ; a þ 1:96 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 À aÞ=R p Þ as the 95% confidence interval for the estimated type I error rate for the nominal level α.