A Novel Approach Integrating Hierarchical Clustering and Weighted Combination for Association Study of Multiple Phenotypes and a Genetic Variant

As a pivotal research tool, genome-wide association study has successfully identified numerous genetic variants underlying distinct diseases. However, these identified genetic variants only explain a small proportion of the phenotypic variation for certain diseases, suggesting that there are still more genetic signals to be detected. One of the reasons may be that one-phenotype one-variant association study is not so efficient in detecting variants of weak effects. Nowadays, it is increasingly worth noting that joint analysis of multiple phenotypes may boost the statistical power to detect pathogenic variants with weak genetic effects on complex diseases, providing more clues for their underlying biology mechanisms. So a Weighted Combination of multiple phenotypes following Hierarchical Clustering method (WCHC) is proposed for simultaneously analyzing multiple phenotypes in association studies. A series of simulations are conducted, and the results show that WCHC is either the most powerful method or comparable with the most powerful competitor in most of the simulation scenarios. Additionally, we evaluated the performance of WCHC in its application to the obesity-related phenotypes from Atherosclerosis Risk in Communities, and several associated variants are reported.

Introduction

Traditionally, Genome-Wide Association Studies studies (GWASs) aim to identify genetic variants associated with certain phenotypes for explaining complex diseases (O’Reilly et al., 2012; Yang and Wang, 2012). In GWASs, multiple related phenotypes of diseases are typically collected for getting better understand complex diseases (Yang Q. et al., 2010). For example, hypertension is directly dependent on the magnitudes of Systolic Blood Pressures (SBP) and Diastolic Blood Pressures (DBP) (Yang and Wang, 2012). The correlation coefficient between SBP and DBP is greater than 0.5 in 95% of patients (Gavish et al., 2008), and researchers could acquire SBP and DBP together for studying hypertension. Similarly, Type 2 Diabetes (T2D) study often gathers relevant risk factors and a number of diabetes-related quantitative phenotypes. Moreover, metabolic syndrome encompasses insulin resistance, obesity, atherosclerotic dyslipidemia, and hypertension; and these factors are interrelated to share potential genetic mediators, pathways, and mechanisms (Huang, 2009). In statistical genetics, jointly analyzing multiple phenotypes can enhance the power of association tests to identify genetic markers associated with one or more phenotypes (Aschard et al., 2014). One of the common approaches for analyzing multiple related phenotypes is to conduct single-phenotype separately and report the results for each phenotype (O’Reilly et al., 2012). However, analysis for one phenotype at a time will be inevitably subject to multiple testing corrections, which leads to a power loss in GWAS (Yang Q. et al., 2010). In recent years, joint analysis of multiple phenotypes has become catching on because of its enhanced statistical power in the detection of genetic variants compared to analysis for each phenotype separately (Yang Q. et al., 2010; Aschard et al., 2014).

Nowadays, jointly analyzing multiple phenotypes has been developed rapidly, which is of two categories: univariate analysis and multivariate analysis. Univariate analysis, as the name suggests, combines various test statistics or p-values of univariate association analysis by means of different strategies. Recently, some approaches of univariate analysis have been proposed for exploring the association between multiple phenotypes and a genetic variant (van der Sluis et al., 2013; Liang et al., 2016; Yang et al., 2016). For example, Kwak et al. (Kwak and Pan, 2016) established an adaptive testing approach, which employs summary statistics from GWASs to evaluate the relationship between multiple phenotypes and a genetic variant. TATES mainly conducts p-values from the association between phenotypes and Single Nucleotide Polymorphisms (SNPs) and concurrently adjusting the correlations among various phenotypes (van der Sluis et al., 2013). Adaptive Fisher’s combination (AFC) (Liang et al., 2016) combines a number of optimal p-values from the traditional GWASs. Compared to multivariate analysis, univariate analysis is generally in a unified framework and tends to ignore the crucial information among multiple phenotypes, which may result in reducing statistical power. In recent years, a series of multivariate analysis approaches including mixed-effects models (Korte et al., 2012; Zhou and Stephens, 2014; Casale et al., 2015), Generalized Estimating Equation (GEE) (Zeger and Liang, 1986; Zhang et al., 2014), and reverse regression methods (O’Reilly et al., 2012; Yan et al., 2013; Wang et al., 2016) have been developed. Mixed-effects models comprise Linear Mixed Effects model (LME) model and Generalized Linear Mixed effects Model (GLMM), where the genetic variants are regarded as the fixed effects and the correlation among phenotypes is considered as random effects. The GEE method collapses the random effects and random residual errors in marginal regression models, which makes it different from LME. The reverse regression methods regard genotypes as the response variable and multiple phenotypes as predictors, such as the proportional odds logistic regression for joint model of multiple phenotypes (MultiPhen) (O’Reilly et al., 2012). Multivariate analysis methods are complicated, and few available software has been developed to implement these methods (Yang and Wang, 2012).

In this article, we establish a novel allele-based approach aiming at detecting association between multiple phenotypes and a genetic variant for improving the power in association studies. We first employ the Hierarchical Clustering based on Different methods for calculating Correlation coefficients (HCDC) (Fu, 2020) to cluster the enrolled phenotypes into several groups. Then, inspired by Weighted Combination of multiple Phenotypes (WCmulP) (Zhu et al., 2018), which provides optimal weights in linear combination, we perform WCmulP in each cluster to generate a novel phenotype by virtual of combining the multiple phenotypes. Subsequently, for every cluster, score test derived from the logistic regression model is constructed to test the association between the genetic variant and the novel phenotype. The permutation procedure is employed to evaluate the p-values of the score test statistics, and their minimum is taken as the test statistic for detecting association between the genetic variant and all phenotypes. Consequently, the Weighted Combination of multiple phenotypes following Hierarchical Clustering method (WCHC) is proposed. Using extensive simulation scenarios, we compare the performance of WCHC with that of six existing methods: O’Brien’s method (O’Brien, 1984), MultiPhen (O’Reilly et al., 2012), MANOVA (Cole et al., 1994), SHet (Zhu et al., 2015), TATES (van der Sluis et al., 2013), and WCmulP (Zhu et al., 2018). The results reveal that WCHC is either the most powerful test or comparable with the most powerful tests among the methods we compared in most of the simulation scenarios. Finally, we evaluate the performance of WCHC approach by utilizing the obesity-related phenotypes from a real dataset, Atherosclerosis Risk in Communities (ARIC) Study from dbGaP, and 11 obesity-associated SNPs are detected.

Materials and Methods

Proposed WCHC

Suppose a sample of N individuals each have M quantitative phenotypes Y₁,Y₂,…,Y_M and genotype G at a genetic variant. It is straightforward to calculate the correlation coefficient between two sets of phenotypes. Based on our previous work (Fu, 2020), the hierarchical clustering is conducted, and finally we have K clusters C₁,C₂,…,C_K. Let M_k denote the number of phenotypes in the kth cluster C_k, k = 1,2,…,K. We take the first cluster C₁ as an example to show the subsequent procedure. Without loss of generality, assume Y₁,Y₂,…,Y_M1 are the M₁ phenotypes in the first cluster. Borrowing the allele-based regression idea (Zhu et al., 2018), we introduce x_2i–1 = x_2i = 1, x_2i–1 = x_2i = 0, and x_2i–1 = 1 and x_2i = 0, if the genotype of the ith individual is AA, aa, and Aa, respectively, i = 1,2,…,N. By analogy, let y_2i,j = y_2i−1,j be the value of the jth phenotype of individual i, i = 1,2,…,N, j = 1,2,…,M₁. Based on ${\{x_l,y_{l,1},y_{l,2},\mathrm{\dots },y_{l,M_1}\}}_{l=1}^{2N}$, we establish the following model:

$$\text{logit}P(x_l=1|y_{l,1},y_{l,2},\mathrm{\dots },y_{l,M_1})={\mathrm{\beta }}_0+{\mathrm{\beta }}_1{y}_{l,1}+{\mathrm{\beta }}_2{y}_{l,2}+\mathrm{\cdots }+{\mathrm{\beta }}_{M_1}{y}_{l,M_1},l=1,2,\mathrm{\dots },2N$$

to test the association between multiple phenotypes Y₁,Y₂,…,Y_M1 and a genetic variant.

Instead of the conventional score test that is vulnerable in the case of big M₁, we adopt the following test statistic (Zhu et al., 2018):

$$T_1=\sum _{l=1}^{2N}(x_l-\overline{x})(y_l-\overline{y}),$$

where $\overline{x}=\frac{1}{2N}{\sum }_{l=1}^{2N}{x}_l$, $\overline{y}=\frac{1}{2N}{\sum }_{l=1}^{2N}{y}_l$, $y_l={\sum }_{j=1}^{M_1}{w}_j{y}_{l,j}$, $w_j=\frac{{\sum }_{l=1}^{2N}(x_l-\overline{x})(y_{l,j}-{\overline{y}}_j)}{{\sum }_{l=1}^{2N}{(y_{l,j}-{\overline{y}}_j)}^2}$, ${\overline{y}}_j=\frac{1}{2N}{\sum }_{l=1}^{2N}{y}_{l,j}$, j = 1,2,…,M₁. Similarly, we have the corresponding test statistics T₂,…,T_K when we study the association of the genetic variant and the multiple phenotypes in the clusters C₂,…,C_K, respectively. Further, let p₁,p₂,…,p_K be the p-values of T₁,T₂,…,T_K, respectively, and we propose our test statistic:

$$T_{WCHC}=\min \{p_1,p_2,\mathrm{\dots },p_K\}.$$

As it is not easy to derive the distribution of the test statistics T₁,T₂,…,T_K under the null hypothesis of no association, the permutation procedure described below is employed to calculate the p-value of T_WCHC.

(1) In each of the B permutations, we random shuffle the genotypes and then get the statistics $T_1^{(b)},T_2^{(b)}\mathrm{\dots },T_K^{(b)}$, b = 0,1,2,…,B. Note that b = 0 is corresponding to the original data (no permutation).

(2) Calculate $p_k^{(b)}$ by:

$p_k^{(b)}=\frac{\mathrm{\#}\{d:T_k^{(d)}>T_k^{(b)}\text{for}d=0,1,\mathrm{\dots },B\}}{B}$, for k = 1,2,…,K.

and then $T_{WCHC}^{(b)}=\min \{p_1^{(b)},p_2^{(b)},\mathrm{\dots },p_K^{(b)}\}$ for b = 0,1,2,…,B;

(3) Then, the p-value of T_WCHC is given by:

$$\frac{\mathrm{\#}\{b:T_{WCHC}^{(b)}<T_{WCHC}^{(0)}\text{for}b=1,2,\mathrm{\dots },B\}}{B}.$$

The hierarchical clustering based on our previous work (Fu, 2020) is as follows: In summary, we can find a partition ψ that partitioned M phenotypes into K disjoint clusters C₁,C₂,…,C_K, where ψ = {C₁,C₂,…,C_K} with ${\bigcup }_{k=1}^K{C}_k=\{1,2,\mathrm{\dots },M\}$ and C_k⋂C_l = ∅(k≠l). Specifically, applying the bottom-up hierarchical clustering approach, we begin with each phenotype as a singleton cluster and then subsequently merge pairs of clusters with the largest similarity until all clusters have been merged into a single cluster that contains all phenotypes. The largest similarity in each iteration is referred as the height of the merged cluster in the dendrogram. A stopping criterion determines the number of clusters, which is similar to an established principle (Bühlmann et al., 2013). Suppose h_b is the largest similarity between two clusters in iteration b (b ≥ 1) or the height of iteration b. We define:

$$\widehat{b}=\arg \underset{b\ge 1}{min}\left(h_{b+1}-h_b\right).$$

Then, the number of clusters identified at the iteration $\widehat{b}$ is chosen to determine the K clusters C₁,C₂,…,C_K. On the calculation of correlation coefficient, the Pearson correlation coefficient, multiple correlation coefficient, and canonical correlation coefficient are respectively employed according to the number of phenotypes in the merged two clusters.

The source code for WCHC method can be found in https://github.com/YQHuFD/WCHC.

Comparison of Methods

For convenience, let 1_n be the all ones vector of length n and 0_n be the all zeroes vector of length n, where n is a positive integer. We first list the following existing methods for power comparison with the proposed WCHC.

OB (O’Brien’s method) (O’Brien, 1984): Using a linear combination of univariate statistics, the OB statistic, ${\mathbf{1}}_M^T{\mathbf{\Sigma }}^{-1}{T}_{uni}$, is developed. It is the most powerful statistic when a class of statistics is a linear combination of T_uni, where T_uni is the vector of univariate statistics and Σ is the variance–covariance matrix of T_uni.

MultiPhen (Joint model of Multiple Phenotypes) (O’Reilly et al., 2012): Modeling the genotype data as ordinal response and phenotypes as predictors, MultiPhen employs likelihood ratio test to evaluate the null hypothesis in the proportional odds logistic regression.

MANOVA (Multivariate ANalysis Of Variance) (Cole et al., 1994): In the standard MANOVA, there are a total of M phenotypes, and the M × M symmetrical background variance–covariance matrix Σ is unconstrained. It has ((M + 1) × M)/2 freely estimated elements in covariances and variances. Standard MANOVA tests the null hypothesis that the M regression coefficients are all zeroes, which is asymptotically equal to the F-test.

SHet (Test for Heterogeneous genetic effects) (Zhu et al., 2015): The test statistic of SHet, S_Het, is based on S_Hom, which is the most powerful statistic when the genetic effects are homogeneous. $S_{Hom}=\frac{{\mathbf{1}}_M^T{(CorrW)}^{-1}{T}_{uni}{({\mathbf{1}}_M^T{(CorrW)}^{-1}{T}_{uni})}^T}{{\mathbf{1}}_M^T{(WCorrW)}^{-1}{\mathbf{1}}_M}$, where Corr is the correlation matrix of T_uni, W is a diagonal matrix of weights for the univariate statistic. S_Het is the maximum of S_Hom’s satisfying various thresholds. Specifically, only the statistics with absolute values greater than the given threshold are employed; Corr and W are partially used corresponding to the selected statistics. The p-value of S_Het could be estimated by simulation.

TATES (Trait-based Association Test that uses Extended Simes procedure) (van der Sluis et al., 2013): TATES combines the p-values of univariate analysis for getting a comprehensive p-value, while correcting the correlation between phenotypes. The TATES p-value is denoted as $\min \left(\frac{M_e{p}_{(j)}}{M_{e(j)}}\right)$, where p_(j) is the j^th (j = 1,…,M) sorted p-value in ascending order; M_e and M_e(j) denote the effective number of independent p-values among all M phenotypes and m specific phenotypes, respectively. The effective numbers can be obtained from the correlation matrix of p-values.

WCmulP (Weighted Combination of multiple Phenotypes) (Zhu et al., 2018): WCmulP can be taken as a component of WCHC. The original phenotypes are not used clustering and directly applied the logistic regression. Then the T statistic is proposed to test the association between the phenotypes and genetic variants. Lastly, the permutation procedure is used to derive the distribution of the test statistic T.

Simulation Studies

Assume that the population is in Hardy–Weinberg equilibrium (HWE), and the genotypes of the genetic variants follow the binomial distribution with parameter 2 and the minor allele frequency (MAF). We set MAF = 0.3 in this simulation study for all scenarios. The multiple phenotypes are generated via the following factor model (van der Sluis et al., 2013):

$$y=\mathrm{\lambda }x+c\mathrm{\gamma }f+\sqrt{1-c^2}\times \mathrm{\epsilon },$$

where y = (y₁,…,y_M)^T is the M phenotypes; x is the genotype; λ = (λ₁,…,λ_M)^T is the vector of effect sizes of the variant on the M phenotypes; f is the vector of factors; f = (f₁,…,f_R)^T∼MVN(0,Σ),Σ = (1−ρ)I + ρA; I is the identity matrix; A is a matrix with elements of 1; R is the number of factors; and ρ is the correlation between factors; γ is an M × R matrix; c is a constant; ε = (ε₁,…,ε_M)^T is a vector of random errors; and ε₁,…,ε_M are mutually independent and follow the standard normal distributions. Consider the following six models with varied numbers of factors.

Model 1: There is only one factor, and the genotype has an influence on all phenotypes with the same effect size. Namely, R = 1, λ=β1_M, and γ=1_M.

Model 2: There are two factors and a genotype has an effect on one factor with the same effect. That is, R = 2, $\mathrm{\lambda }={({\mathbf{0}}_{M\mathrm{/}2}^T,\mathrm{\beta }{\mathbf{1}}_{M\mathrm{/}2}^T)}^T$, and γ = bdiag(1_M/2,1_M/2), which represents the block diagonal matrix of 1_M/2 and 1_M/2.

Model 3: There are two factors, and a genotype has an effect on the second factor with different sizes. That is, R = 2, $\mathrm{\lambda }={({\mathbf{0}}_{M\mathrm{/}2}^T,\frac{\mathrm{\beta }}{M+1}{[1:M\mathrm{/}2]}^T+\mathrm{\beta }{\mathbf{1}}_{M\mathrm{/}2}^T)}^T$ and γ = bdiag(1_M/2,1_M/2), where [1:M/2] represents the vector of components 1,2,…,M/2.

Model 4: There are four factors, and a genotype has an impact on the last factor with the same size. That is, R = 4, $\mathrm{\lambda }={({\mathbf{0}}_{3M\mathrm{/}4}^T,\mathrm{\beta }{\mathbf{1}}_{M\mathrm{/}4}^T)}^T$, and γ = bdiag(1_M/4,1_M/4,1_M/4,1_M/4).

Model 5: There are four factors, and a genotype has an effect on the last factor with different sizes. Namely, R = 4, $\mathrm{\lambda }={({\mathbf{0}}_{3M\mathrm{/}4}^T,\frac{\mathrm{\beta }}{M+1}{[1:M\mathrm{/}4]}^T+\mathrm{\beta }{\mathbf{1}}_{M\mathrm{/}4}^T)}^T$, γ = bdiag(1_M/4,1_M/4,1_M/4,1_M/4).

Model 6: There are four factors, and a genotype has an impact on the last two factors with different effect directions. That is, R = 4, $\mathrm{\lambda }={({\mathbf{0}}_{M\mathrm{/}2}^T,-\frac{\mathrm{\beta }}{M+1}{[1:M\mathrm{/}4]}^T-\mathrm{\beta }{\mathbf{1}}_{M\mathrm{/}4}^T,\mathrm{\beta }{\mathbf{1}}_{M\mathrm{/}4}^T)}^T$, γ = bdiag(1_M/4,1_M/4,1_M/4,1_M/4).

For these six models, the within-factor correlation is c² and the between-factor correlation is ρc². For estimating type I error rates and powers, we fix N = 1,000 unrelated subjects, the number of phenotypes M = 16, 32. By means of setting β = 0, we generate all phenotypes that is independent of genotypes to evaluate the type I error rates of all methods, including OB, MultiPhen, MANOVA, SHet, TATES, WCmulP, and WCHC. The corresponding Q–Q plot of type I error rates is shown in Supplementary Figures 1–6. Importantly, to evaluate powers, we not only vary the values of β (while within-factor correlation c² = 0.5 and between-factor correlation ρc² = 0.1) but also change the values of within-factor correlation c² = 0.3, 0.5, 0.7, and 0.9 (while between-factor correlation ρc² = 0.1).

The calculation of heritability is as follows: the heritability of genotypes to the j-th phenotype is given by $h^2(y_j)=\frac{var(x){\mathrm{\lambda }}_j^2}{var(x){\mathrm{\lambda }}_j^2+1}\approx var(x){\mathrm{\lambda }}_j^2.$ The heritability of genotypes to the total M phenotypes is given by $h^2={\sum }_{j=1}^M{h}^2(y_j)\approx var(x){\sum }_{j=1}^M{\mathrm{\lambda }}_j^2.$ Then given the parameters λ, M, and MAF, we can calculate h² for the different models.

Simulation Results

We set different nominal significance levels, various numbers of phenotypes, and distinct number of factors to estimate the type I error rates of WCHC and other six methods. For each simulation scenario, the p-values of WCHC, WCmulP, and SHet are evaluated by 2,000 permutations; and the p-values of MANOVA, MultiPhen, TATES, and OB are evaluated by their asymptotic distributions. The type I error rates of the seven methods are estimated using 2,000 replicated samples. For 2,000 replicated samples, the 95% confidence intervals (CIs) for type I error rates of nominal levels 0.01 and 0.05 are about (0.0056, 0.0144) and (0.0404, 0.0596), respectively. The evaluated type I error rates of WCHC and other six methods are presented in Table 1 (M = 16) and Table 2 (M = 32). It is observed from these two tables that most of the type I error rates of WCHC are within 95% CIs, which shows the validity of the developed WCHC. Meanwhile, the type I error rates of WCmulP, SHet, MANOVA, MultiPhen, TATES, and OB are not obviously deviated from the nominal levels. See more information in Q–Q plots (Supplementary Figures 1–6).

TABLE 1

Table 1. Type I error rates of the seven methods in three simulation settings.

TABLE 2

Table 2. Type I error rates of the seven methods in three simulation settings.

In order to compare powers of these seven methods, we plot power against the genetic effect β (in Figures 1, 2) and the within-factor correlation c² (in Figures 3, 4). Note in the calculation of power, the p-values of WCHC, WCmulP, and SHet are evaluated by 1,000 permutations; the powers of the seven methods are estimated based on 1,000 replicated samples at a significance level of 0.05. The following observations can be drawn from the simulation.

FIGURE 1

Figure 1. Power comparisons of the seven methods as a function of β in the six models. Sample size is N = 1,000, the number of phenotypes is M = 16, c² = 0.5, ρc² = 0.1, and MAF = 0.3. The power of all the seven methods is estimated by 1000 replicated samples at a significance level of 0.05.

FIGURE 2

Figure 2. Power comparisons of the seven methods as a function of β in the six models. Sample size is N = 1,000, the number of phenotypes is M = 32, c² = 0.5, ρc² = 0.1, and MAF = 0.3. The power of all the seven methods is estimated by 1000 replicated samples at a significance level of 0.05.

FIGURE 3

Figure 3. Power comparisons of the seven methods as a function of c² in the six models. Sample size is N = 1,000, the number of phenotypes is M = 16, ρc² = 0.1 and MAF = 0.3. β = 0.09 for model 1 and 2; β = 0.08 for model 3; β = 0.1 for model 4 and 5; β = 0.07 for model 6. The power of all the seven methods is estimated by 1000 replicated samples at a significance level of 0.05.

FIGURE 4

Figure 4. Power comparisons of the seven methods as a function of c² in the six models. Sample size is N = 1,000, the number of phenotypes is M = 32, ρc² = 0.1, and MAF = 0.3. β = 0.1 for model 1 and 4–6; β = 0.09 for model 2; β = 0.08 for model 3. The power of all the seven methods is estimated by 1000 replicated samples at a significance level of 0.05.

(1) As expected, in each model, the powers of all seven methods increase as the genetic effect β increases (see Figures 1, 2). (2) Except in models 1 and 6, WCHC is the most powerful test in all the methods under most of the simulation scenarios (see Figures 1–4). (3) As the number of phenotypes increases from M = 16 to M = 32, WCHC exhibits more obvious advantages over other methods except in Model 1 and 6 (see Figures 1, 2). (4) No matter changes of genetic effects β or variations of correlation coefficients between different phenotypes, MANOVA and MultiPhen have the similar performance in all six models. (5) Generally, in each model, the power of all methods decreases with the increase of correlation coefficients of within factors between phenotypes. (6) OB is the most powerful test when the genetic effects are homogeneous (Model 1). However, OB’s power decreases when there exist opposite directions (Model 6) or when the genetic variant has an influence on a small proportion of phenotypes (Model 5). (7). In general, WCHC, WCmulP, and TATES are more powerful than SHet, OB, MANOVA, and MultiPhen when the genetic variant affects a portion of phenotypes (Models 2–6). (8). WCHC shows obvious advantages over other methods when the genetic variant only affects part of the phenotypes with the same direction. One possible reason is that in the models of generating data, the genetic variant has effects of the same directions on some phenotypes and has no effect on the remaining ones. The hierarchical clustering is capable of grouping similar phenotypes together, so as to reduce the dimensions of association test for improving the power to detect the associated phenotypes.

Overall, from all the power simulation results, we could draw that our proposed WCHC has advantages over other methods in most scenarios, and especially in some scenes, the ascendancy is obvious. In other scenarios, WCHC is comparable with the most powerful test.

Real Data Analysis

We applied our proposed method WCHC to the real data analysis from ARIC study (see more details in The ARIC Investigators, 1989). In brief, sponsored by the National Heart, Lung, and Blood Institute (NHLBI), ARIC is a prospective cohort study of atherosclerosis risk in community. It records the changes of the incidence of atherosclerosis-related diseases and cardiovascular risk factors in distinct races, regions, genders, and time, aiming at investigating the etiology and natural process of atherosclerosis (Morrison et al., 2013). We obtain the genotyped and clinical phenotypic data in ARIC from dbGaP server of the United States National Center for Biotechnology Information (accession number: phs000090.v4.p1).

To evaluate the performance of WCHC in real data, we use the seven methods to analyze obesity-related phenotypes in ARIC. We selected nine continuous traits with regard to obesity including weight, body mass index (BMI), average skinfold thickness of the triceps brachii, mean subscapular skinfold thickness, waist, hip girth, waist-to-hip ratio, calf girth, and wrist breadth and three covariates including age, gender, and race. The specific description of these variables is listed in Table 3, and the correlation structure of obesity-related phenotypes is given in Supplementary Figure 7. A set of 12,701 subjects across 272,027 SNPs were left for subsequent analysis after excluding subjects with missing data in any of the 12 variables as well as the genetic variants with missing rate greater than 0.2 or HWE < 10^–4. Every phenotype is adjusted for those three covariates using linear regression model.

TABLE 3

Table 3. Summary statistics of obesity-related indexes and covariates in ARIC.

Based on these adjusted phenotypes related to obesity, we employ WCHC and other six methods to detect associated SNPs. Two groups are obtained after clustering the nine phenotypes in the real data analyses by the hierarchical clustering in WCHC. One of the clusters only includes wrist breadth, while the other encompasses the remaining phenotypes. Because of multiple testing correction, we adopt the significance threshold of 1 × 10^–7, not the traditional genome-wide significance threshold of 5 × 10^–8. There are totally 11 SNPs that are significant for at least one method (Table 4). Previous studies (Frayling et al., 2007; Heard-Costa et al., 2009; Lindgren et al., 2009; Meyre et al., 2009; Thorleifsson et al., 2009; Willer et al., 2009; Heid et al., 2010; Speliotes et al., 2010; Bradfield et al., 2012; Wen et al., 2012; Berndt et al., 2013; Monda et al., 2013; Locke et al., 2015; Shungin et al., 2015) have reported that FTO leads to obesity through population studies and experimental researches elaborating relevant mechanisms. Among the 11 identified SNPs, rs9939609 and rs8050136 are involved in FTO. Additionally, rs7968682 is reported to be associated with height (Yang T. L. et al., 2010; Takeshita et al., 2011). Few other SNPs have been assessed to explore the association with obesity or obesity-related phenotypes. From Table 4, we can see that both WCHC and MANOVA identified six SNPs; TATES identified five SNPs; both WCmulP and SHet identified four SNPs; MultiPhen identified three SNPs; and OB only identified one SNP, which may be due to that the true genetic effects of most of SNPs are heterogeneous for all phenotypes. In summary, the number of SNPs identified by WCHC is comparable with the largest number of SNPs identified by other tests. These real data analysis results are consistent with our simulation results.

TABLE 4

Table 4. Significant SNPs and the corresponding p-values in the analysis of ARIC.

Characteristics of the Significant Variants

Table 5 shows the annotations of the identified SNPs based on the Ensemble website¹ and SCAN website². From Table 5, we can see that the significant SNPs are located in intergenic or intron region, and most of them have been reported to be associated with BMI, height, or T2D. Generally, they have been reported in GWAS. We could also explore the expression of genes associated with the significant SNPs, although they are located in intergenic or intron region. Therefore, we make full use of Qtlizer³, eQTLGen⁴, and PsychENCODE⁵, which are the largest integrating various tissues, blood, and brain expression Quantitative Trait Locus (eQTL) samples, respectively. Instead of restricting analysis to the SNPs in Table 5, we considered using a larger list of SNPs with proxy variants, which are in Linkage Disequilibrium (LD) with the SNPs in Table 5 (r² ≥ 0.8) via Qtlizer website. We restricted the eQTL association criteria with False Discovery rate (FDR) < 0.05. The results of eQTLs in Qtlizer, eQTLGen, and PsychENCODE are displayed in Supplementary Data Sheets.

TABLE 5

Table 5. Characteristics of the significant SNPs.

In order to further study the biology function of the genetic variants, we performed enrichment analysis on genes associated with these 11 SNPs in Table 5 and the proxy variants (see qtlizer.results in Supplementary Data Sheet) in the three websites/consortiums (Qtlizer, eQTLGen, and PsychENCODE). After summarizing all genes in the three tables (see qtlizer.results, eQTLGen.results, and PsychENCODE.results in Supplementary Data Sheet), we got all the genes associated with the eQTLs (see Gene sheet in Supplementary Data Sheet). A total of 76 genes were obtained to do the gene set analysis by virtue of different biological databases for investigating biological processes, cell components, molecular functions, metabolic pathways, phenotypes with relevant diseases, and protein interactions. The results of enrichment analysis and protein–protein interaction (PPI) are given in Figures 5, 6. According to the Gene Ontology (GO) enrichment analysis chart in Figure 5, GO items mainly focus on the cellular response to hydrogen and regulation of lipid kinase activity, which may be parts of the metabolic process. Moreover, the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway in Figure 6 presents that the enriched genes possess taurine and hypotaurine metabolism and endocrine resistance, which indicates that the obesity-related variants detected by WCHC and other methods in ARIC might be involved in the metabolic pathways, which releases the signal that our results are in a certain degree of credibility. Subsequently, we draw a PPI network diagram through the STRING⁶ to reveal that most of the proteins (nodes) encoded by the genes have certain interactions (edges), which suggests the proteins related the expression of genes might interact with each other for controlling a variety of biological phenomena including endocrine development, cellular response to hydrogen, and metabolic processes.

FIGURE 5

Figure 5. GO enrichment analysis of significant SNPs probability regulating associated genes expression. (A) Red, blue, and green bars indicate biology progress, cellular components, and molecular function categories, respectively. The numbers above the bar charts indicate the number of genes in each of the biological categories; (B) Bar charts of GO enrichment analysis; (C) Volcano plot of GO enrichment analysis. For more knowledge about GO enrichment, please check the website http://geneontology.org/docs/go-enrichment-analysis/.

FIGURE 6

Figure 6. KEGG enrichment analysis and PPI network diagram of significant SNPs probability regulating associated genes expression. (A) Bar chats of KEGG enrichment analysis; (B) Volcano plot of KEGG enrichment analysis; (C) PPI interaction network diagram, data are from https://www.string-db.org/.

Overall, our results showed that WCHC and other six methods could identify significant genetic variants for obesity phenotypes in real data analysis from ARIC. More importantly, functional annotations of genetic variants and enrichment analysis support that the variants are closely related to biological functions and metabolic pathways of obesity.

Discussion

In this article, we proposed WCHC to perform multivariate analysis of multiple phenotypes in association studies due to the following reasons. (1) Multiple correlated phenotypes are usually measured in complex disease for genetic association studies. Compared to univariate analysis, multivariate analysis considers multidimensional structure information. It indicates certain variance–covariance is included in multiple phenotypes. (2) Association analysis of multiple phenotypes separately cannot present genetic interactions between phenotypes. More and more evidence reveals that joint analysis of multiple related phenotypes, which considers the interactions between phenotypes comprehensively, can boost the power of detecting genetic variants associated with complex diseases. No matter whether the effects of genetic variants on phenotypes are consistent or not, WCHC provides a relatively simple way to incorporate the correlations between phenotypes into analysis. (3) Actually, we are not sure which phenotype or linear combination of phenotypes is more likely to elucidate the genetic structure of complex diseases. WCHC adopts clustering approach and linear combination of multiple phenotypes to account for the complex genetic information, which not only takes the similarity between phenotypes into consideration but also considers the heterogeneity, so it is helpful to explore the genetic mechanism of diseases.

Our results manifested that WCHC has correct type I error rates and is either the most powerful test or comparable with the most powerful tests among the seven methods we adopted. None of the other methods observes consistently good performances under the simulation scenarios. OB is the most powerful test when the genetic effects are homogeneous, while it loses power dramatically when genetic effects are heterogeneous, especially if there exists opposite directions of genetic effects. In most simulation scenarios, SHet, MANOVA, MultiPhen, and TATES have similar powers, and they are less powerful than WCHC, and WCHC is more powerful when the genetic variant influences a part of phenotypes. However, WCmulP is less powerful in this scenario. Furthermore, in real data analysis, WCHC and MANOVA identified the largest number of significant SNPs (six SNPs). Therefore, the real data analysis results demonstrate that WCHC has excellent performance in detecting SNPs associated with complex disease with multiple related phenotypes such as obesity. As for the methods giving such different results when applied to the real ARIC data, we think that the parametric information of real data is unknown for us. Therefore, we may try various methods to analyze the real data for getting reliable results as much as possible. In our opinion, no method can guarantee 100% accuracy. We can only be cautious to say that the significant loci are more likely to be true signals, but further verification is still needed.

In the context of association studies, population stratification (PS) refers to allele frequency difference between populations uncorrelated to the outcome of interest, but due to systematic ancestry differences. PS may cause confounding effects seriously if not adjusted properly (Knowler et al., 1988; Lander and Schork, 1994). Methods such as principal component analysis (PCA) (Zhu et al., 2002; Chen et al., 2003; Zhang et al., 2003; Price et al., 2006; Bauchet et al., 2007), linear mixed model (LMM) (Kang et al., 2010; Zhang et al., 2010; Hoffman, 2013), multidimensional scaling (MDS) (Li and Yu, 2008), robust PCA based on resampling by half means (RPCA-RHM) (Liu et al., 2013), and robust PCA based on the projection pursuit (RPCA-PP) (Liu et al., 2013) can be used to adjust for PS. We propose to apply PCA to control for PS when samples from different populations are involved.

In real data analysis, as the number of phenotypes elevates, the chance of missing at least one subject increases exponentially, especially in epidemiological and clinical research (Ali et al., 2011; Dahl et al., 2016). We removed 412 subjects with missing either phenotypes or covariates from 13,113 observations. It is worth noting that the sample mean substitution (Ali et al., 2011; van der Sluis et al., 2013) is a simple, unconditional method that does not depend on other variables, which is a common strategy replacing the missing values with plausible values for the variable with missing values. However, it may contribute to biased estimates where data are not missing completely at random (Ali et al., 2011). Additionally, imputation is a more complicated approach that fills in missing values with estimated values via model-based methods or conditional imputation, comprising multiple imputation (MI), multivariate normal imputation (MVNI), and fully conditional specification (FCS) (Raghunathan et al., 2001; Buuren et al., 2006; De Silva et al., 2017).

One weakness of WCHC is that the test statistic does not have an asymptotic distribution and its p-value needs to be calculated by permutation procedure, which is time-consuming as compared with approaches whose test statistics have asymptotic distributions. To conduct GWAS, a small number of permutations (e.g., 1,000) can be used to select genetic variants that reveal evidence of association, and then a large number of permutations are employed to estimate the selected significant genetic variants. We adopted this strategy to analyze the real dataset. Consequently, it seems to be efficient, and the bioinformatics analysis of significant variants supports our results. In conclusion, in the field of genotype–phenotype association studies, WCHC is an effective method for association analysis of multiple phenotypes, which considers both the correlations and differences among the multiple phenotypes. WCHC provides a convenient approach of association analysis for researchers to discover potential genes causing complex diseases, which does not need to assume the genetic model, and there is no limit to the number of phenotypes. Because the genetic structure of phenotypes is usually unknown, WCHC provides a convenient statistical method for the application of massive multi-phenotypic data in the future.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: phs000090.v4.p1. The source code for WCHC method can be found in https://github.com/YQHuFD/WCHC.

Ethics Statement

The studies involving human participants were reviewed and approved by the ARIC. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

LF and Y-QH: study concept and design and drafting of the manuscript. LF, Y-QH, and YW: acquisition of data. LF: methodology and interpretation of data. LF, YW, TL, and Y-QH: critical revision of the manuscript for important intellectual content. All authors have read and approved the final version of manuscript. We thank all reviewers and editors for their valuable suggestions on revision.

Funding

This study was supported by grants to Y-QH from the National Natural Science Foundation of China (grant nos. 11971117 and 11571082).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.654804/full#supplementary-material

Footnotes

1. ^ https://asia.ensembl.org
2. ^ http://scandb.org
3. ^ http://genehopper.de/qtlizer
4. ^ https://www.eqtlgen.org/
5. ^ http://resource.psychencode.org/
6. ^ https://www.string-db.org/

References

Ali, A. M., Dawson, S. J., Blows, F. M., Provenzano, E., Ellis, I. O., Baglietto, L., et al. (2011). Comparison of methods for handling missing data on immunohistochemical markers in survival analysis of breast cancer. Br. J. Cancer 104, 693–699. doi: 10.1038/sj.bjc.6606078

PubMed Abstract | CrossRef Full Text | Google Scholar

Aschard, H., Vilhjálmsson, B. J., Greliche, N., Morange, P.-E., Trégouët, D.-A., and Kraft, P. (2014). Maximizing the power of principal-component analysis of correlated phenotypes in genome-wide association studies. Am. J. Hum. Genet. 94, 662–676. doi: 10.1016/j.ajhg.2014.03.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bauchet, M., McEvoy, B., Pearson, L. N., Quillen, E. E., Sarkisian, T., Hovhannesyan, K., et al. (2007). Measuring European population stratification with microarray genotype data. Am. J. Hum. Genet. 80, 948–956. doi: 10.1086/513477

PubMed Abstract | CrossRef Full Text | Google Scholar

Berndt, S. I., Gustafsson, S., Mägi, R., Ganna, A., Wheeler, E., Feitosa, M. F., et al. (2013). Genome-wide meta-analysis identifies 11 new loci for anthropometric traits and provides insights into genetic architecture. Nat. Genet. 45, 501–512. doi: 10.1038/ng.2606

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradfield, J. P., Taal, H. R., Timpson, N. J., Scherag, A., Lecoeur, C., Warrington, N. M., et al. (2012). A genome-wide association meta-analysis identifies new childhood obesity loci. Nat. Genet. 44, 526–531. doi: 10.1038/ng.2247

PubMed Abstract | CrossRef Full Text | Google Scholar

Bühlmann, P., Rutimann, P., van de Geer, S., and Zhang, C. H. (2013). Correlated variables in regression: clustering and sparse estimation. J. Stat. Plan Infer. 143, 1835–1858. doi: 10.1016/j.jspi.2013.05.019

CrossRef Full Text | Google Scholar

Buuren, S., Brand, J., Groothuis-Oudshoorn, C., and Rubin, D. (2006). Fully conditional specification in multivariate imputation. J. Stat. Comput. Simul. 76, 1049–1064. doi: 10.1080/10629360600810434

CrossRef Full Text | Google Scholar

Casale, F. P., Rakitsch, B., Lippert, C., and Stegle, O. (2015). Efficient set tests for the genetic analysis of correlated traits. Nat. Methods 12, 755–758. doi: 10.1038/nmeth.3439

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H. S., Zhu, X., Zhao, H., and Zhang, S. (2003). Qualitative semi-parametric test for genetic associations in case-control designs under structured populations. Ann. Hum. Genet. 67, 250–264. doi: 10.1046/j.1469-1809.2003.00036.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cole, D. A., Maxwell, S. E., Arvey, R., and Salas, E. (1994). How the power of MANOVA can both increase and decrease as a function of the intercorrelations among the dependent variables. Psychol. Bull. 115, 465–474.

Google Scholar

Dahl, A., Iotchkova, V., Baud, A., Johansson, Å, Gyllensten, U., Soranzo, N., et al. (2016). A multiple-phenotype imputation method for genetic studies. Nat. Genet. 48, 466–472. doi: 10.1038/ng.3513

PubMed Abstract | CrossRef Full Text | Google Scholar

De Silva, A. P., Moreno-Betancur, M., De Livera, A. M., Lee, K. J., and Simpson, J. A. (2017). A comparison of multiple imputation methods for handling missing values in longitudinal data in the presence of a time-varying covariate with a non-linear association with time: a simulation study. BMC Med. Res. Methodol 17:114. doi: 10.1186/s12874-017-0372-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Frayling, T. M., Timpson, N. J., Weedon, M. N., Zeggini, E., Freathy, R. M., Lindgren, C. M., et al. (2007). A common variant in the FTO gene is associated with body mass index and predisposes to childhood and adult obesity. Science 316, 889–894. doi: 10.1126/science.1141634

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, L. (2020). Multi-Phenotype Association studies and Mendelian Randomization for Pleiotropic Genetic Variants and their Applications. Ph.D. thesis. Shanghai: Fudan Univesity.

Google Scholar

Gavish, B., Ben-Dov, I. Z., and Bursztyn, M. (2008). Linear relationship between systolic and diastolic blood pressure monitored over 24 h: assessment and correlates. J. Hypertens. 26, 199–209. doi: 10.1097/HJH.0b013e3282f25b5a

PubMed Abstract | CrossRef Full Text | Google Scholar

Heard-Costa, N. L., Zillikens, M. C., Monda, K. L., Johansson, A., Harris, T. B., Fu, M., et al. (2009). NRXN3 is a novel locus for waist circumference: a genome-wide association study from the charge consortium. PLoS Genet. 5:e1000539. doi: 10.1371/journal.pgen.1000539

PubMed Abstract | CrossRef Full Text | Google Scholar

Heid, I. M., Jackson, A. U., Randall, J. C., Winkler, T. W., Qi, L., Steinthorsdottir, V., et al. (2010). Meta-analysis identifies 13 new loci associated with waist-hip ratio and reveals sexual dimorphism in the genetic basis of fat distribution. Nat. Genet. 42, 949–960. doi: 10.1038/ng.685

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffman, G. E. (2013). Correcting for population structure and kinship using the linear mixed model: theory and extensions. PLoS One 8:e75707. doi: 10.1371/journal.pone.0075707

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, P. L. (2009). A comprehensive definition for metabolic syndrome. Dis. Model. Mech. 2, 231–237. doi: 10.1242/dmm.001180

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, H. M., Sul, J. H., Service, S. K., Zaitlen, N. A., Kong, S.-Y., Freimer, N. B., et al. (2010). Variance component model to account for sample structure in genome-wide association studies. Nat. Genet. 42, 348–354. doi: 10.1038/ng.548

PubMed Abstract | CrossRef Full Text | Google Scholar

Knowler, W. C., Williams, R. C., Pettitt, D. J., and Steinberg, A. G. (1988). Gm3;5,13,14 and type 2 diabetes mellitus: an association in American Indians with genetic admixture. Am. J. Hum. Genet. 43, 520–526.

Google Scholar

Korte, A., Vilhjálmsson, B. J., Segura, V., Platt, A., Long, Q., and Nordborg, M. (2012). A mixed-model approach for genome-wide association studies of correlated traits in structured populations. Nat. Genet. 44, 1066–1071. doi: 10.1038/ng.2376

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwak, I. Y., and Pan, W. (2016). Adaptive gene- and pathway-trait association testing with GWAS summary statistics. Bioinformatics 32, 1178–1184. doi: 10.1093/bioinformatics/btv719

PubMed Abstract | CrossRef Full Text | Google Scholar

Lander, E. S., and Schork, N. J. (1994). Genetic dissection of complex traits. Science 265, 2037–2048.

Google Scholar

Li, Q., and Yu, K. (2008). Improved correction for population stratification in genome-wide association studies by identifying hidden population structures. Genet. Epidemiol. 32, 215–226. doi: 10.1002/gepi.20296

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, X., Wang, Z., Sha, Q., and Zhang, S. (2016). An adaptive Fisher’s combination method for joint analysis of multiple phenotypes in association studies. Sci. Rep. 6:34323. doi: 10.1038/srep34323

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindgren, C. M., Heid, I. M., Randall, J. C., Lamina, C., Steinthorsdottir, V., Qi, L., et al. (2009). Genome-wide association scan meta-analysis identifies three loci influencing adiposity and fat distribution. PLoS Genet. 5:e1000508. doi: 10.1371/journal.pgen.1000508

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Zhang, D., Liu, H., and Arendt, C. (2013). Robust methods for population stratification in genome wide association studies. BMC Bioinformatics 14:132. doi: 10.1186/1471-2105-14-132

PubMed Abstract | CrossRef Full Text | Google Scholar

Locke, A. E., Kahali, B., Berndt, S. I., Justice, A. E., Pers, T. H., Day, F. R., et al. (2015). Genetic studies of body mass index yield new insights for obesity biology. Nature 518, 197–206. doi: 10.1038/nature14177

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyre, D., Delplanque, J., Chèvre, J.-C., Lecoeur, C., Lobbens, S., Gallina, S., et al. (2009). Genome-wide association study for early-onset and morbid adult obesity identifies three new risk loci in European populations. Nat. Genet. 41, 157–159. doi: 10.1038/ng.301

PubMed Abstract | CrossRef Full Text | Google Scholar

Monda, K. L., Chen, G. K., Taylor, K. C., Palmer, C., Edwards, T. L., Lange, L. A., et al. (2013). A meta-analysis identifies new loci associated with body mass index in individuals of African ancestry. Nat. Genet. 45, 690–696. doi: 10.1038/ng.2608

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrison, A. C., Voorman, A., Johnson, A. D., Liu, X., Yu, J., Li, A., et al. (2013). Whole-genome sequence-based analysis of high-density lipoprotein cholesterol. Nat. Genet. 45, 899–901. doi: 10.1038/ng.2671

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Brien, P. C. (1984). Procedures for comparing samples with multiple endpoints. Biometrics 40, 1079–1087.

Google Scholar

O’Reilly, P. F., Hoggart, C. J., Pomyen, Y., Calboli, F. C., Elliott, P., Jarvelin, M. R., et al. (2012). Multiphen: joint model of multiple phenotypes can increase discovery in GWAS. PLoS One 7:e34861. doi: 10.1371/journal.pone.0034861

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A., and Reich, D. (2006). Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909. doi: 10.1038/ng1847

PubMed Abstract | CrossRef Full Text | Google Scholar

Raghunathan, T., Lepkowski, J., Hoewyk, J., and Solenberger, P. (2001). A multivariate technique for multiply imputing missing values using a sequence of regression models. Surv. Methodol 27, 85–95.

Google Scholar

Shungin, D., Winkler, T. W., Croteau-Chonka, D. C., Ferreira, T., Locke, A. E., Mägi, R., et al. (2015). New genetic loci link adipose and insulin biology to body fat distribution. Nature 518, 187–196. doi: 10.1038/nature14132

PubMed Abstract | CrossRef Full Text | Google Scholar

Speliotes, E. K., Willer, C. J., Berndt, S. I., Monda, K. L., Thorleifsson, G., Jackson, A. U., et al. (2010). Association analyses of 249,796 individuals reveal 18 new loci associated with body mass index. Nat. Genet. 42, 937–948. doi: 10.1038/ng.686

PubMed Abstract | CrossRef Full Text | Google Scholar

Takeshita, H., Fujihara, J., Soejima, M., Koda, Y., Kimura-Kataoka, K., Ono, R., et al. (2011). Confirmation that SNPs in the high mobility group-a2 gene (HMGA2) are associated with adult height in the Japanese population; wide-ranging population survey of height-related SNPs in HMGA2. Electrophoresis 32, 1844–1851. doi: 10.1002/elps.201100128

PubMed Abstract | CrossRef Full Text | Google Scholar

The ARIC Investigators (1989). The atherosclerosis risk in communities (ARIC) study: design and objectives. Am. J. Epidemiol. 129, 687–702.

Google Scholar

Thorleifsson, G., Walters, G. B., Gudbjartsson, D. F., Steinthorsdottir, V., Sulem, P., Helgadottir, A., et al. (2009). Genome-wide association yields new sequence variants at seven loci that associate with measures of obesity. Nat. Genet. 41, 18–24. doi: 10.1038/ng.274

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Sluis, S., Posthuma, D., and Dolan, C. V. (2013). TATES: efficient multivariate genotype-phenotype analysis for genome-wide association studies. PLoS Genet. 9:e1003235. doi: 10.1371/journal.pgen.1003235

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Sha, Q., and Zhang, S. (2016). Joint analysis of multiple traits using “optimal” maximum heritability test. PLoS One 11:e0150975. doi: 10.1371/journal.pone.0150975

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, W., Cho, Y.-S., Zheng, W., Dorajoo, R., Kato, N., Qi, L., et al. (2012). Meta-analysis identifies common variants associated with body mass index in east Asians. Nat. Genet. 44, 307–311. doi: 10.1038/ng.1087

PubMed Abstract | CrossRef Full Text | Google Scholar

Willer, C. J., Speliotes, E. K., Loos, R. J. F., Li, S., Lindgren, C. M., Heid, I. M., et al. (2009). Six new loci associated with body mass index highlight a neuronal influence on body weight regulation. Nat. Genet. 41, 25–34. doi: 10.1038/ng.287

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, T., Li, Q., Li, Y., Li, Z., and Zheng, G. (2013). Genetic association with multiple traits in the presence of population stratification. Genet. Epidemiol. 37, 571–580. doi: 10.1002/gepi.21738

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. J., Li, J., Williams, L. K., and Buu, A. (2016). An efficient genome-wide association test for multivariate phenotypes based on the Fisher combination function. BMC Bioinformatics 17:19. doi: 10.1186/s12859-015-0868-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Q., and Wang, Y. (2012). Methods for analyzing multivariate phenotypes in genetic association studies. J. Probab. Stat. 2012:652569. doi: 10.1155/2012/652569

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Q., Wu, H., Guo, C. Y., and Fox, C. S. (2010). Analyze multivariate phenotypes in genetic association studies by combining univariate association tests. Genet. Epidemiol. 34, 444–454. doi: 10.1002/gepi.20497

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, T. L., Guo, Y., Zhang, L. S., Tian, Q., Yan, H., Guo, Y. F., et al. (2010). HMGA2 is confirmed to be associated with human adult height. Ann. Hum. Genet. 74, 11–16. doi: 10.1111/j.1469-1809.2009.00555.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeger, S. L., and Liang, K. Y. (1986). Longitudinal data analysis for discrete and continuous outcomes. Biometrics 42, 121–130.

Google Scholar

Zhang, S., Zhu, X., and Zhao, H. (2003). On a semiparametric test to detect associations between quantitative traits and candidate genes using unrelated individuals. Genet. Epidemiol. 24, 44–56. doi: 10.1002/gepi.10196

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Xu, Z., Shen, X., and Pan, W. (2014). Testing for association with multiple traits in generalized estimation equations, with application to neuroimaging data. Neuroimage 96, 309–325. doi: 10.1016/j.neuroimage.2014.03.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Ersoz, E., Lai, C.-Q., Todhunter, R. J., Tiwari, H. K., Gore, M. A., et al. (2010). Mixed linear model approach adapted for genome-wide association studies. Nat. Genet. 42, 355–360. doi: 10.1038/ng.546

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., and Stephens, M. (2014). Efficient multivariate linear mixed model algorithms for genome-wide association studies. Nat. Methods 11, 407–409. doi: 10.1038/nmeth.2848

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, H., Zhang, S., and Sha, Q. (2018). A novel method to test associations between a weighted combination of phenotypes and genetic variants. PLoS One 13:e0190788. doi: 10.1371/journal.pone.0190788

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Feng, T., Tayo, B. O., Liang, J., Young, J. H., Franceschini, N., 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. doi: 10.1016/j.ajhg.2014.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Zhang, S., Zhao, H., and Cooper, R. S. (2002). Association mapping, using a mixture model for complex traits. Genet. Epidemiol. 23, 181–196. doi: 10.1002/gepi.210

PubMed Abstract | CrossRef Full Text | Google Scholar