Associating Multivariate Traits with Genetic Variants Using Collapsing and Kernel Methods with Pedigree- or Population-Based Studies

Chien, Li-Chu

doi:https://doi.org/10.1155/2021/8812282

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Materials and Methods Results Conclusion Data Availability Conflicts of Interest Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2021 | Article ID 8812282 | https://doi.org/10.1155/2021/8812282

Associating Multivariate Traits with Genetic Variants Using Collapsing and Kernel Methods with Pedigree- or Population-Based Studies

Li-Chu Chien¹

Academic Editor: Giuseppe Pontrelli

Received06 Sept 2020

Revised02 Jan 2021

Accepted08 Jan 2021

Published10 Feb 2021

Abstract

In genetic association analysis, several relevant phenotypes or multivariate traits with different types of components are usually collected to study complex or multifactorial diseases. Over the past few years, jointly testing for association between multivariate traits and multiple genetic variants has become more popular because it can increase statistical power to identify causal genes in pedigree- or population-based studies. However, most of the existing methods mainly focus on testing genetic variants associated with multiple continuous phenotypes. In this investigation, we develop a framework for identifying the pleiotropic effects of genetic variants on multivariate traits by using collapsing and kernel methods with pedigree- or population-structured data. The proposed framework is applicable to the burden test, the kernel test, and the omnibus test for autosomes and the X chromosome. The proposed multivariate trait association methods can accommodate continuous phenotypes or binary phenotypes and further can adjust for covariates. Simulation studies show that the performance of our methods is satisfactory with respect to the empirical type I error rates and power rates in comparison with the existing methods.

1. Introduction

Genome-wide association studies (GWAS) intend to find genetic variants such as single nucleotide polymorphisms (SNPs) associated with common traits or with complex diseases [1, 2]. Association studies, where the correlation relationship between a genetic variant and a trait is evaluated, are helpful for mapping genes influencing complex diseases [3]. In the study of complex diseases, data on several correlated phenotypes or a multivariate phenotype with several components are often collected to get a better understanding of the disease [1, 3, 4]. Multivariate correlated traits are influenced through multiple variants simultaneously. Therefore, by a suitable joint or multivariate analysis framework of multivariate traits, we can not only gain more statistical power to identify pleiotropic effects of genetic variants on multivariate traits [3, 5–12] but also can further understand the genetic architecture of the disease of interest [5, 13]. Thus, recently, the joint analysis of multivariate traits has become popular because it can increase statistical power over analyzing only one trait at a time [1, 4].

Several statistical methods have been developed to identify the association between multivariate traits and a genetic variant [1, 5]. Current multivariate methods can be classified into three groups [1, 2, 5]: regression methods [14–16], variable reduction methods [11, 13, 17, 18], and combining analysis [9, 19–23]. However, many of the existing methods for multivariate association analysis cannot be straightaway extended to rare variant analyses, due to their enormous numbers causing the problems of multiple comparison or multiple testing and their low minor allele frequencies [2, 5, 24]. Moreover, sparsity of data could lead to problems on estimating regression parameters and fitting regression models [2]. Hence, it is necessary for proposing statistical methods for identifying the association between multivariate traits and multiple genetic variants (common and/or rare variants) [5]. In recent years, various statistical techniques have been proposed for this purpose in GWAS [8, 17, 25–27]. Furthermore, several approaches have been extendedly developed for the investigation of rare variants associated with multivariate traits [2, 28–38].

Although these new developments keep many benefits, existing methods have some potential limitations [39]. Most current methods are constructed under some specific assumptions about the effects of genetic variants on multivariate traits [39]. These current approaches suffer a severe loss in power once the model assumptions are violated [26, 39].

In this investigation, we develop the statistical methods for identifying pleiotropic effects of genetic variants on multivariate traits using collapsing and kernel methods with pedigree- or population-structured data. The proposed multivariate trait association method is able to handle binary phenotypes or continuous phenotypes and further can adjust for covariates. Moreover, the proposed multivariate trait association method not only can leverage the dependence on the phenotypes but also can account for the sample relatedness in the pedigree-based or population-based structured data.

The rest of the article is organized as follows. In the materials and methods section, we construct the multivariate effect model using the joint GEE model formulation (JGEE) [40]. We apply the JGEE to pedigree- or population-structured data and introduce a retrospective framework for analyzing multivariate traits in genetic association studies. The proposed framework is applicable to the burden test, the kernel test, and the omnibus test for autosomes and the X chromosome. In the simulation studies, we examine the finite sample size performance of the proposed multivariate association methods and evaluate the comparison results with the existing method, Multi-SKAT [39]. Concluding remarks and future possibilities for continuity are given in the conclusion section and the limitation section.

2. Materials and Methods

2.1. Notations

To describe the proposed multivariate trait association method based on the pedigree- or population-based structured data, we suppose that there are independent pedigrees and each pedigree has subjects. We assume that the subjects have been sequenced in a genetic region of interest (e.g., a gene) that contains variants. Let be the phenotype vector for the ^th phenotype of the ^th pedigree. Let be the response vector for the phenotypes that we are interested in. Let be the vector for the ^th covariate of the ^th pedigree. Let be the covariate matrix for the nongenetic covariates that we want to adjust for. Let be the vector of regression coefficients of the nongenetic covariates with the element being the effect of the ^th covariate on the ^th trait. Let be the genetic matrix for genetic variants in a target region of interest where is the vector for a genetic variant (, 1, or 2 for 0, 1, or 2 copies of the minor allele, respectively). Let be the vector of regression coefficients of the genetic variants with the element being the effect of the ^th genetic variant on the ^th trait.

2.2. Multitrait Regression-Based Tests for Pedigree Data

We let be the covariate matrix and be the genotype matrix for the ^th pedigree where is an identity matrix of dimension and stands for the Kronecker product. According to the generalized linear model [41], we assume that the marginal density of is with two moments, and , where is a scale parameter. Let be the vector with the components and be the vector with the components for the ^th trait of the ^th pedigree.

Based on the joint GEE model formulation [40], we construct the multivariate linear model for describing the association relationship between correlated traits and genetic variants, which is given as follows: where is the inverse function of and is a response-specific link function [40], is the vector of the expected mean of the multivariate traits , is the vector of regression coefficients of the nongenetic covariates for the correlated traits, and is the vector of regression coefficients of the genetic variants for the correlated traits.

Let and be the within-in cluster correlation matrix and the multivariate-response cluster correlation matrix, which depend on a vector of parameters and , respectively. The working (or approximate) covariance matrix of is given by [40]. where is a block diagonal matrix with the components being the diagonal matrices. According to equation (1), under the null hypothesis of no association between genotypes and phenotypes, we propose the multivariate association methods including the homogeneous kernel statistic (HoK), the heterogeneous kernel statistic (HeK), and burden test (BT). Moreover, we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO) by combining the HoK with the BT and by combining the HeK with the BT, respectively.

2.2.1. Kernel Statistic

We let be a correlation matrix of genotype scores with element for markers and . Let denote the minor allele frequency (MAF) of the ^th marker. Let be the vector of the standard residuals with components where is the inverse matrix of . Here, and are the estimates of and . Here and henceforth, all estimates are calculated based on the null hypothesis of the genetic effects equal to zero. All unknown parameters and the working within-in and multivariate-response cluster correlation matrices are estimated by the R package JGEE [42].

(1) Homogeneous Kernel Statistic. We suppose that is a marker-specific weight of the ^th variant and assume that the genetic effects on the different phenotypes are homogeneous (i.e., . Based on the JGEE model with the genotype as random variables considered, we propose the homogeneous quadratic (kernel) association statistic (HoK) as follows: where , , _, is the estimate of , and is the estimate of that is a diagonal matrix for the ^th phenotype of the ^th pedigree. The null distribution of asymptotically follows a mixture chi-square distribution , where s are independent random variables following a chi-square distribution with one degree of freedom and are nonzero eigenvalues of the null covariate matrix of where and is a matrix of genetic correlations for all individuals in the ^th pedigree, which has the same definition given by Schaid et al. [43] and can be calculated by the R package kinship2 [44]. When the genetic relationship between subjects and in the ^th pedigree is unknown, the elements of the genetic correlation can be estimated through genomic data [43, 45], and its estimate is given by [43]

(2) Heterogeneous Kernel Statistic. We assume that the genetic effects on the different phenotypes are heterogeneous (i.e., . The heterogeneous quadratic (kernel) association statistic (HeK) is defined by where and is a marker-specific weight of the ^th variant of the ^th trait. The null distribution of asymptotically follows a mixture chi-square distribution , where s are independent random variables following a chi-square distribution with one degree of freedom, and are nonzero eigenvalues of the null covariate matrix of , where .

Theoretical values of and are approximately calculated by Kuonen’s saddlepoint method [46] and obtained by the R package pchisqsum. A theory for the derivation of the HoK test and the HeK test is shown in Appendix S1.

2.2.2. Burden Test

We let be a weighted average of genotype scores for the ^th pedigree. On the basis of the HoK test and the HeK test in equations (3) and (5) with the same marker-specific weight of the ^th variant for each trait (i.e., ), we propose the burden test (BT) as follows: where the null covariance matrix of is given by

Then,

The null distribution of asymptotically follows a chi-square distribution with one degree of freedom.

2.2.3. Omnibus Test

Let , , and denote the values obtained by the HoK, HeK, and BT statistics. Based on the idea of the value combination method through the Cauchy distribution [47–49], we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO).

(1) Homogeneous Omnibus Test. Combining the with the , we construct the homogeneous omnibus test (HoO) as follows: where stands for the inverse cumulative distribution function of the standard Cauchy distribution.

(2) Heterogeneous Omnibus Test. Combining the with the , we construct the heterogeneous omnibus test (HeO) as follows:

The null distributions of the test and the test asymptotically follow a standard Cauchy distribution [47–49]. The values of the test and the test are calculated by the R package RNOmni [50].

The kernel statistic, the burden test, and the omnibus test are also applicable to the X chromosome. Additional technical information for extensions to the X chromosome is shown in Appendix S2.

3. Simulation Studies

We conduct the numerical simulation studies to assess the finite sample performance of the proposed methods and evaluate the comparison results with two existing methods, the minimum value SKAT statistic (mPK), and the minimum value burden statistic (mPB) [39]. The two existing methods are implemented by the R package Multi-SKAT [39]. Based on the similar simulation set-up as those usually considered from existing genetic association tests [39, 43, 51], we investigate the effect of the proposed methods, HoK, HeK, BT, HoO, and HeO, for identifying genetic variants that are associated with multiple traits. We simultaneously generate 10,000 European-like (EUR) and 10,000 admixed African American-like (AA) haplotypes of length 200 kb using a calibrated human demographic model through the COSI software [51, 52]. A 3 kb region is randomly selected in our numerical simulations. We generate a total of 10,000 databases for each simulation scenario in our studies.

3.1. Type I Error Rate and Power Simulations

In the heterogeneous population with nuclear family data considered, continuous and binary phenotypes for trait for individual in the ^th family are generated from the multivariate linear model in equation (1) with and . More precisely, continuous and binary phenotypes are generated by the following linear and logit models, respectively: where , , , , , and . Here, the elements of the covariance matrix is a vector of all ones. The elements of are independently generated with an equal probability of being 0 or 1. The elements of are generated from a multivariate normal distribution with a mean of 0.5 and a covariance matrix with diagonal entries of 1 and all off-diagonal entries of 0.1. The regression coefficients of the covariate matrix for the ^th correlated trait are given by and , respectively, for continuous traits and binary traits for .

For continuous traits, the error terms in equation (11) follow a multivariate normal distribution having a mean of zero, a within-in cluster correlation matrix (i.e., ) with diagonal entries of 1 and all off-diagonal entries of 0.2 and a subject-across-response correlation matrix (i.e., ) with diagonal entries of 0.3 and all off-diagonal entries of 0.1. Similarly, binary traits in equation (12) are generated with the same within-in cluster correlation matrix (i.e., ) and the same subject-across-response correlation matrix (i.e., ) as the continuous traits in equation (11). These correlated phenotypes are generated by the R package BinNor [53].

For type I error simulations, the regression coefficients of genetic variants, , in equations (11) and (12) are equal to zero under the null hypothesis. For power simulations, under the alternative hypothesis, we simulate that 35% of low variants with the are causal. For each setting, either all causal SNPs have a positive effect, or 80% of causal SNPs are positive, and 20% of causal SNPs are negative. The regression coefficients of genetic variants, , are set by or corresponding to the risk or protective variant [51]. Under the assumption that the genetic effects on the two different phenotypes are heterogeneous (i.e., , the genetic effects for the first traits are set as described above, while the genetic effects for the second traits are set by zero. On the other hand, under the assumption that the genetic effects on the two different phenotypes are homogeneous (i.e., , the genetic effects and for the first and second traits have the same settings as described above.

We simulate 1,400 nuclear families with 800 nuclear families from the European samples and 600 nuclear families from African-American samples. The marker-specific weight for variant is considered as the beta density function with shape parameters and [51]. To study the effect of the marker-specific weight of variant on the phenotypes, we consider the unweighted marker-specific weight with and the weighted marker-specific weight with [51]. The empirical type I error rates based on fifty thousand replicates and the empirical power rates based on two thousand replicates are reported for all simulation results. The “exchangeable” and “unstructured” structures are considered for the working within-cluster and multivariate-response correlation matrices for the proposed methods, HoK, HeK, and BT, respectively.

4. Results

4.1. Empirical Type I Error Rates

Table 1 reports the results of a simulation comparison on empirical type I error rates when the phenotypes are considered to be continuous. Table 1 displays that the proposed methods, HoK, HoO, HeK, HeO, and BT, well control the empirical type I error rates regardless of the weight of the marker-specific weight. Similarly, the existing methods, mPK and mPB, have good performance on controlling the empirical type I error rates. Our simulation results show that the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, reasonably control the empirical type I error rates for autosome analyses with continuous traits. The seven competing approaches display similar performance in terms of the empirical type I error rates for the X chromosome analyses with continuous traits (Appendix S3: Table S1).

Table 2 reports the empirical type I error rates based on the proposed methods, HoK, HeK, BT, HoO, and HeO, for the binary data. The two existing methods, mPK and mPB, aren’t included for comparison. This reason is that implementing the two existing methods, mPK and mPB, via the R package Multi-SKAT [39], the MPMM (multiple phenotype mixed model) function in the R package PHENIX [54–56] is a necessary tool for this process. However, the MPMM function is suitable for the continuous phenotypes [56] or is suitable for the binary phenotypes with the condition that the number of cases is sufficiently large [39]. In other words, in some sense, the two existing methods, mPK and mPB, are limited to continuous phenotypes [39].

Table 2 shows that the proposed methods appropriately control the type I error rates when the marker-specific weight is considered for or for variant for binary traits. On the other hand, the empirical type I error rates of the proposed methods for X chromosome analyses with binary traits are depicted in Table S2 in Appendix S3. These empirical type I error rates show similar results as that for autosome analyses.

In summary, our simulation results show that the proposed multivariate trait association methods, HoK, HoO, HeK, HeO, and BT, have reasonable control of type I error rates for continuous traits or binary traits whether the marker is X chromosomal or autosomal. On the other hand, the existing methods, mPK and mPB, yield well-controlled type I error rates for the autosome analyses or the X chromosome analyses with continuous traits (Table 1 or Table S1), regardless of the weight of the marker-specific weight.

4.2. Empirical Power

Figure 1 exhibits the comparison results of the empirical power rates for the autosome analyses with continuous traits, when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be exchangeable. As expected, the empirical power rates of the seven competing methods with a weighted marker-specific weight of are higher than that with an unweighted marker-specific weight of . The heterogeneous kernel statistic (HeK) has slightly greater empirical power rates than other methods, when the genetic effects on the different phenotypes are heterogeneous (i.e., ), and causal SNPs have positive effects or negative effects on phenotypes. On the other hand, the existing method, mPB, has bigger empirical power rates, when the genetic effects on the different phenotypes are heterogeneous (i.e., ), and all causal SNPs have a positive association on phenotypes. Moreover, the empirical power rates of the homogeneous omnibus test (HoO) are larger than that of the other six competing methods, when the genetic effects on the different phenotypes are homogeneous (i.e., ). Evidently, the seven competing methods have their respective advantages in identifying the association between genetic effects and multiple continuous traits for autosome analyses.

Similar empirical power rates are obtained from the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, considered to be unstructured. Hence, these empirical power rates are not shown in order to save space. On the other hand, the seven competing approaches display a similar performance in testing for the X chromosome analyses with continuous traits (Appendix S3: Figure S1).

Figure 2 exhibits the comparison results of empirical power rates for the autosome analyses with binary traits when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be exchangeable. As a similar reason for investigating the empirical type I error rates with binary traits, the two existing methods, mPK and mPB, aren’t included for power comparison.

Figure 2 shows that the heterogeneous kernel statistic (HeK) and the heterogeneous omnibus test (HeO) outperform over other methods in terms of the empirical power rates, when the genetic effects on the different phenotypes are heterogeneous (i.e., ). On the other hand, the empirical power rates of the homogeneous omnibus test (HoO) are bigger than that of the other competing methods, when the genetic effects on the different phenotypes are homogeneous (i.e., ). As expected, in general, the heterogeneous kernel statistic (HeK) is more powerful than the homogeneous kernel statistic (HoK), when the genetic effects on the different phenotypes are heterogeneous (i.e., ). On the other hand, the homogeneous kernel statistic (HoK) is more powerful than the heterogeneous kernel statistic (HeK), when the genetic effects on the different phenotypes are homogeneous (i.e., ). In a word, the proposed methods, HoK, HoO, HeK, HeO, and BT, have their respective merits in examining the association between genetic effects and multiple binary traits for autosome analyses.

Similarly, when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be unstructured, the empirical power rates have similar results and thus they are omitted. On the other hand, the empirical power rates of the proposed methods for X chromosome analyses with binary traits are presented in Figure S2 in Appendix S3. These empirical power rates show similar results as that discussed in Figure 2.

In summary, the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, have their respective merits in diagnosing whether genetic effects are associated with multiple continuous traits for autosome analyses or the X chromosome analyses. Similarly, the proposed methods, HoK, HoO, HeK, HeO, and BT, have their respective advantages in examining whether there are associations between genetic effects and multiple binary traits for autosome analyses or the X chromosome analyses.

To furthermore examine the performance of the proposed methods, additional simulation studies for continuous traits and binary traits are presented in Appendix S4 and Appendix S5 with higher correlations of phenotypes and higher dimensions of phenotypes considered, respectively. In general, these competing methods based on higher correlations of phenotypes or higher dimensions of phenotypes can provide a bigger empirical power rate for the analysis of continuous traits or binary traits. However, we note that these competing methods based on higher correlations of phenotypes or higher dimensions of phenotypes more easily have empirical type I error rate inflation at a smaller nominal level, especially for binary data analysis (Appendix S5: Tables S5-S6 and Appendix S6: Table S7), in comparison with these methods based on lower correlations of phenotypes or lower dimensions of phenotypes. A detailed discussion of these additional simulation results is given in Appendixes S4 and S5.

However, we note that the proposed methods have a high computational cost, especially for binary data. Under our simulation setting and framework, we carry out a single simulated data set by using a computer based on one CPU core at 2.1 GHz. The average computational times of the homogeneous and heterogeneous tests with a weighted marker-specific weight under the alternative hypothesis for continuous data are 0.83 and 0.91 minutes, respectively, while that for binary data are 4.77 and 4.80 minutes, respectively. Therefore, in the current version, such a framework algorithm implementation is unsatisfactory for analyzing a large-scale high-dimensional data set in practice.

5. Conclusion

In this investigation, we develop a retrospective framework for identifying the pleiotropic effects of genetic variants on multivariate traits by using collapsing and kernel methods with pedigree- or population-structured data. The proposed framework, corresponding to the burden test, the kernel test, and the omnibus test, provides a sound basis for genetic association analyses for autosomes and the X chromosome. The proposed multivariate trait association methods based on the JGEE model can flexibly accommodate continuous phenotypes or binary phenotypes and further can adjust for covariates.

One critical advantage of the proposed methods is that the homogeneous kernel statistic (HoK), the heterogeneous kernel statistic (HeK), and the burden test (BT) retain all of the benefits of the retrospective tests proposed by Schaid et al. [43] who treated the genotype data as random variables by conditioning the phenotypes as constants. On the other hand, the homogeneous omnibus test (HoO) and the heterogeneous kernel statistic (HeO) keep the advantages of the Cauchy combination tests proposed by Liu and Xie [48] who showed that the Cauchy combination tests are robust to model misspecification and robustly protect the type I error rates [49].

Another important benefit of the proposed method is that the HoK test, the HeK test, and the BT test keep the benefits of the JGEE model that validly account for complex correlations between subjects within the cluster (within-cluster correlations) and between different phenotypes from the same subjects (multivariate-response correlations). Moreover, the proposed test statistics, HoK, HeK, and BT, based on the JGEE model can efficaciously account for covariate adjustment whether the phenotypes are continuous or binary.

Our simulation studies show that an unweighted marker-specific weight and an exchangeable structure of the working within-cluster and multivariate-response correlations are recommended for the practical data analysis if the data cannot sufficiently provide valid information for estimating the structures of the working within-cluster and multivariate-response correlations before the start of the data analysis. Moreover, the homogeneous kernel statistic (HoK) is more robust than the heterogeneous statistic (HeK) in controlling the empirical type I errors, because the null distribution of the HeK statistic asymptotically follows a mixture chi-square distribution with a larger degree of freedom, in comparison with the null distribution of the HoK statistic. However, the HeK statistic is more powerful than the HoK statistic when the genetic effects on the different phenotypes are heterogeneous.

On the other hand, our simulation results show that for the autosome analyses or the X chromosome analyses with continuous traits, the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, show good performance with well-controlled type I errors, while the seven competing methods have their respective merits for identifying the association between the genetic effects and multiple continuous traits. In addition, our simulation results show that for the autosome analyses or the X chromosome analyses with binary traits, the proposed methods, HoK, HoO, HeK, HeO, and BT, can control empirical type I errors with lower correlations of phenotypes or with lower dimensions of phenotypes (Table 2 and Table S2), while these proposed methods have their respective advantages for identifying the genetic variants associated with multiple binary traits. However, we observe that the proposed methods, HoK, HoO, HeK, HeO, and BT, with higher correlations of phenotypes or with higher dimensions of phenotypes, more easily have the infection of empirical type I errors at a smaller nominal level (Appendix S5: Tables S5-S6 and Appendix S6: Table S7), although these method under such situations have higher empirical power rates.

6. Limitation

The proposed multivariate trait association methods have their limitations. First, these proposed methods cannot simultaneously include the continuous traits and binary traits in analysis. Thus, future studies are needed to extend the idea of the proposed multivariate trait association methods for simultaneously considering continuous traits and binary traits in analysis. Second, the multivariate trait association methods, based on higher correlations of phenotypes or higher dimensions of phenotypes, easily suffer from the problem of the inflated type I errors, especially when the binary traits are considered (Appendix S5: Tables S5-S6 and Appendix S6: Table S7). Although the JGEE model provides an efficient algorithm for estimating the structure of the working within-cluster and multivariate-response correlations, a large-scale pedigree study always suffers from a more complex and high-dimensional structure of the within-cluster and multivariate-response correlations in pedigree database analysis. Hence, in the future, a more effective algorithm for estimating the complicated and high-dimensional (or higher correlational) structure of the working within-cluster and multivariate-response correlations is necessary to be proposed, especially when the analysis focuses on the binary traits. Third, in comparison with the null distribution of the homogeneous kernel statistic, the null distribution of the heterogeneous kernel statistic follows a larger degree of freedom test, which easily causes such a heterogeneous test to suffer from the problem of the type I error inflation. Therefore, overcoming the problem of the type I error inflation from the heterogeneous test is an essential part of the future work. Fourth, the proposed methods, which have a high computational cost especially for binary data, are inappropriate for analyzing large-scale high-dimensional data in practice. Thus, a more effective algorithm for reducing computational cost is needed to be proposed in further research. Moreover, the software of the proposed methods is computationally inconvenient and particularly inadequate for the mass GWAS data in practice. Therefore, the software of the proposed methods, which is convenient to be used, is a further work in the future. Fifth, our current work focuses mainly on the low- and common-frequency variants. Extension of the proposed methods to the rare variants deserves further works.

Data Availability

The data supporting the findings of this study are available within the article and its supplementary materials.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank the editor and the referees for their constructive comments, which significantly improve the presentation of the article. This work is supported by grant MOST 108-2118-M-037-001-MY2 of Ministry of Science and Technology, Taiwan, R.O.C.

Supplementary Materials

Appendix S1: the null distribution of the kernel statistic. Appendix S2: extension to the X chromosome. Appendix S3: simulation results based on the X chromosome. Appendix S4: additional simulation studies for continuous traits. Appendix S5: additional simulation studies for binary traits. Appendix S6: limitation (Supplementary Materials)

References

H. Zhu, S. Zhang, and Q. Sha, “A novel method to test associations between a weighted combination of phenotypes and genetic variants,” PLoS One, vol. 13, no. 1, article e0190788, 2018.
View at: Publisher Site | Google Scholar
S. Lee, S. Won, Y. J. Kim et al., “Rare variant association test with multiple phenotypes,” Genetic Epidemiology, vol. 41, no. 3, pp. 198–209, 2017.
View at: Publisher Site | Google Scholar
Q. Yang and Y. Wang, “Methods for analyzing multivariate phenotypes in genetic association studies,” Journal of probability and statistics, vol. 2012, Article ID 652569, 13 pages, 2012.
View at: Publisher Site | Google Scholar
X. Liang, Q. Sha, and S. Zhang, “Joint analysis of multiple phenotypes in association studies using allele-based clustering approach for non-normal distributions,” Annals of Human Genetics, vol. 82, no. 6, pp. 389–395, 2018.
View at: Publisher Site | Google Scholar
Z. Wang, Q. Sha, S. Fang, K. Zhang, and S. Zhang, “Testing an optimally weighted combination of common and/or rare variants with multiple traits,” PLoS One, vol. 13, no. 7, article e0201186, 2018.
View at: Publisher Site | Google Scholar
N. Solovieff, C. Cotsapas, P. H. Lee, S. M. Purcell, and J. W. Smoller, “Pleiotropy in complex traits: challenges and strategies,” Nature Reviews Genetics, vol. 14, no. 7, pp. 483–495, 2013.
View at: Publisher Site | Google Scholar
M. Stephens, “A unified framework for association analysis with multiple related phenotypes,” PLoS One, vol. 8, no. 7, article e65245, 2013.
View at: Publisher Site | Google Scholar
X. Zhou and M. Stephens, “Efficient multivariate linear mixed model algorithms for genome-wide association studies,” Nature Methods, vol. 11, no. 4, pp. 407–409, 2014.
View at: Publisher Site | Google Scholar
X. Liang, Z. Wang, Q. Sha, and S. Zhang, “An adaptive Fisher’s combination method for joint analysis of multiple phenotypes in association studies,” Scientific Reports, vol. 6, no. 1, article 34323, 2016.
View at: Publisher Site | Google Scholar
Z. Wang, X. Wang, Q. Sha, and S. Zhang, “Joint analysis of multiple traits in rare variant association studies,” Annals of Human Genetics, vol. 80, no. 3, pp. 162–171, 2016.
View at: Publisher Site | Google Scholar
Z. Wang, Q. Sha, and Z. Zhang, “Joint analysis of multiple traits using "optimal" maximum heritability test,” PLoS One, vol. 11, no. 3, article e0150975, 2016.
View at: Publisher Site | Google Scholar
H. Zhu, S. Zhang, and Q. Sha, “Power comparisons of methods for joint association analysis of multiple phenotypes,” Human heredity, vol. 80, no. 3, pp. 144–152, 2015.
View at: Google Scholar
H. Aschard, B. Vilhjálmsson, N. Greliche, P. Morange, D. Trégouët, and P. Kraft, “Maximizing the power of principal-component analysis of correlated phenotypes in genome-wide association studies,” The American Journal of Human Genetics, vol. 94, no. 5, pp. 662–676, 2014.
View at: Publisher Site | Google Scholar
A. Korte, B. J. Vilhjálmsson, V. Segura, A. Platt, Q. Long, and M. Nordborg, “A mixed-model approach for genome-wide association studies of correlated traits in structured populations,” Nature Genetics, vol. 44, no. 9, pp. 1066–1071, 2012.
View at: Publisher Site | Google Scholar
P. F. O’Reilly, C. J. Hoggart, Y. Pomyen et al., “MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS,” PLoS One, vol. 7, no. 5, article e34861, 2012.
View at: Publisher Site | Google Scholar
Y. Zhang, Z. Xu, X. Shen, W. Pan, and Alzheimer's Disease Neuroimaging Initiative, “Testing for association with multiple traits in generalized estimation equations, with application to neuroimaging data,” NeuroImage, vol. 96, pp. 309–325, 2014.
View at: Publisher Site | Google Scholar
M. A. R. Ferreira and S. M. Purcell, “A multivariate test of association,” Bioinformatics, vol. 25, no. 1, pp. 132-133, 2009.
View at: Publisher Site | Google Scholar
L. Klei, D. Luca, B. Devlin, and K. Roeder, “Pleiotropy and principal components of heritability combine to increase power for association analysis,” Genetic Epidemiology, vol. 32, no. 1, pp. 9–19, 2008.
View at: Publisher Site | Google Scholar
P. C. O’Brien, “Procedures for comparing samples with multiple endpoints,” Biometrics, vol. 40, no. 4, pp. 1079–1087, 1984.
View at: Publisher Site | Google Scholar
Q. Yang, H. Wu, C.-Y. Guo, and C. S. Fox, “Analyze multivariate phenotypes in genetic association studies by combining univariate association tests,” Genetic Epidemiology, vol. 34, no. 5, pp. 444–454, 2010.
View at: Publisher Site | Google Scholar
S. van der Sluis, D. Posthuma, and C. V. Dolan, “TATES: efficient multivariate genotype-phenotype analysis for genome-wide association studies,” PLoS Genet, vol. 9, no. 1, article e1003235, 2013.
View at: Publisher Site | Google Scholar
J. Kim, Y. Bai, and W. Pan, “An adaptive association test for multiple phenotypes with GWAS summary statistics,” Genetic Epidemiology, vol. 39, no. 8, pp. 651–663, 2015.
View at: Publisher Site | Google Scholar
X. Zhu, T. Feng, B. O. Tayo et al., “Meta-analysis of correlated traits via summary statistics from GWASs with an application in hypertension,” The American Journal of Human Genetics, vol. 96, no. 1, pp. 21–36, 2015.
View at: Publisher Site | Google Scholar
B. Li and S. M. Leal, “Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data,” The American Journal of Human Genetics, vol. 83, no. 3, pp. 311–321, 2008.
View at: Publisher Site | Google Scholar
J. Huang, A. Johnson, and C. O'Donnell, “PRIMe: a method for characterization and evaluation of pleiotropic regions from multiple genome-wide association studies,” Bioinformatics, vol. 27, no. 9, pp. 1201–1206, 2011.
View at: Publisher Site | Google Scholar
D. Ray, J. Pankow, and S. Basu, “USAT: a unified score-based association test for multiple phenotype-genotype analysis,” Genetic Epidemiology, vol. 40, no. 1, pp. 20–34, 2016.
View at: Publisher Site | Google Scholar
J. Ried, A. Döring, K. Oexle et al., “PSEA: phenotype set enrichment analysis-a new method for analysis of multiple phenotypes,” Genetic Epidemiology, vol. 36, no. 3, pp. 244–252, 2012.
View at: Publisher Site | Google Scholar
Q. Yan, D. Weeks, J. Celedón et al., “Associating multivariate quantitative phenotypes with genetic variants in family samples with a novel kernel machine regression method,” Genetics, vol. 201, no. 4, pp. 1329–1339, 2015.
View at: Publisher Site | Google Scholar
X. Zhan, N. Zhao, A. Plantinga et al., “Powerful genetic association analysis for common or rare variants with high-dimensional structured traits,” Genetics, vol. 206, no. 4, pp. 1779–1790, 2017.
View at: Publisher Site | Google Scholar
A. Maity, P. Sullivan, and J. Tzeng, “Multivariate phenotype association analysis by marker-set kernel machine regression,” Genetic Epidemiology, vol. 36, no. 7, pp. 686–695, 2012.
View at: Publisher Site | Google Scholar
UK10K Consortium, J. Sun, K. Oualkacha et al., “A method for analyzing multiple continuous phenotypes in rare variant association studies allowing for flexible correlations in variant effects,” European Journal of Human Genetics, vol. 24, no. 9, pp. 1344–1351, 2016.
View at: Publisher Site | Google Scholar
Y. Wang, A. Liu, J. L. Mills et al., “Pleiotropy analysis of quantitative traits at gene level by multivariate functional linear models,” Genetic Epidemiology, vol. 39, no. 4, pp. 259–275, 2015.
View at: Publisher Site | Google Scholar
F. P. Casale, B. Rakitsch, C. Lippert, and O. Stegle, “Efficient set tests for the genetic analysis of correlated traits,” Nature Methods, vol. 12, no. 8, pp. 755–758, 2015.
View at: Publisher Site | Google Scholar
K. A. Broadaway, D. J. Cutler, R. Duncan et al., “A statistical approach for testing cross-phenotype effects of rare variants,” The American Journal of Human Genetics, vol. 98, no. 3, pp. 525–540, 2016.
View at: Publisher Site | Google Scholar
A. Cichonska, J. Rousu, P. Marttinen et al., “metaCCA: summary statistics-based multivariate meta-analysis of genome-wide association studies using canonical correlation analysis,” Bioinformatics, vol. 32, no. 13, pp. 1981–1989, 2016.
View at: Publisher Site | Google Scholar
J. Lin, R. Tabassum, S. Ripatti, and M. Pirinen, “MetaPhat: detecting and decomposing multivariate associations from univariate genome-wide association statistics,” Frontiers in Genetics, vol. 11, 2020.
View at: Publisher Site | Google Scholar
S. Bhattacharjee, P. Rajaraman, K. B. Jacobs et al., “A subset-based approach improves power and interpretation for the combined analysis of genetic association studies of heterogeneous traits,” The American Journal of Human Genetics, vol. 90, no. 5, pp. 821–835, 2012.
View at: Publisher Site | Google Scholar
23andMe Research Team, Social Science Genetic Association Consortium, P. Turley et al., “Multi-trait analysis of genome-wide association summary statistics using MTAG,” Nature Genetics, vol. 50, no. 2, pp. 229–237, 2018.
View at: Publisher Site | Google Scholar
D. Dutta, L. Scott, M. Boehnke, and S. Lee, “Multi-SKAT: general framework to test for rare-variant association with multiple phenotypes,” Genetic Epidemiology, vol. 43, no. 1, pp. 4–23, 2019.
View at: Publisher Site | Google Scholar
G. Inan and R. Yucel, “Joint GEEs for multivariate correlated data with incomplete binary outcomes,” Journal of Applied Statistics, vol. 44, no. 11, pp. 1920–1937, 2017.
View at: Publisher Site | Google Scholar
K. Y. Liang and S. L. Zeger, “Longitudinal data analysis using generalized linear models,” Biometrika, vol. 73, no. 1, pp. 13–22, 1986.
View at: Publisher Site | Google Scholar
G. Inan, “JGEE: joint generalized estimating equation solver,” 2015, R package version 1.1.
View at: Google Scholar
D. J. Schaid, S. K. McDonnell, J. P. Sinnwell, and S. N. Thibodeau, “Multiple genetic variant association testing by collapsing and kernel methods with pedigree or population structured data,” Genetic Epidemiology, vol. 37, no. 5, pp. 409–418, 2013.
View at: Publisher Site | Google Scholar
J. Sinnwell, T. Therneau, D. Schaid, E. Atkinson, and C. Mester, “kinship2,” 2020, R package version 1.8.5.
View at: Google Scholar
T. Thornton and M. S. McPeek, “ROADTRIPS: case-control association testing with partially or completely unknown population and pedigree structure,” The American Journal of Human Genetics, vol. 86, no. 2, pp. 172–184, 2010.
View at: Publisher Site | Google Scholar
H. Chen, J. B. Meigs, and J. Dupuis, “Sequence kernel association test for quantitative traits in family samples,” Genetic Epidemiology, vol. 37, no. 2, pp. 196–204, 2013.
View at: Publisher Site | Google Scholar
Y. Liu, S. Chen, Z. Li, A. C. Morrison, E. Boerwinkle, and X. Lin, “ACAT: a fast and powerful p value combination method for rare-variant analysis in sequencing studies,” The American Journal of Human Genetics, vol. 104, no. 3, pp. 410–421, 2019.
View at: Publisher Site | Google Scholar
Y. Liu and J. Xie, “Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures,” Journal of the American Statistical Association, vol. 115, no. 529, pp. 393–402, 2020.
View at: Publisher Site | Google Scholar
Z. R. McCaw, J. M. Lane, R. Saxena, S. Redline, and X. Lin, “Operating characteristics of the rank-based inverse normal transformation for quantitative trait analysis in genome-wide association studies,” Biometrics, vol. 76, no. 4, pp. 1262–1272, 2020.
View at: Publisher Site | Google Scholar
Z. R. McCaw, “Rank normal transformation omnibus test,” 2019, R package version 0.7.1.
View at: Google Scholar
S. Lee, T. M. Teslovich, M. Boehnke, and X. Lin, “General framework for meta-analysis of rare variants in sequencing association studies,” The American Journal of Human Genetics, vol. 93, no. 1, pp. 42–53, 2013.
View at: Publisher Site | Google Scholar
S. F. Schaffner, C. Foo, S. Gabriel, D. Reich, M. J. Daly, and D. Altshuler, “Calibrating a coalescent simulation of human genome sequence variation,” Genome Research, vol. 15, no. 11, pp. 1576–1583, 2005.
View at: Publisher Site | Google Scholar
A. Amatya, H. Demirtas, and R. Gao, “Simultaneous generation of multivariate binary and normal variates,” 2020, R package version 2.3.2.
View at: Google Scholar
R. Torices and A. J. Muñoz-Pajares, “Phenotypic integration index,” 2017, R package version 1.3.1.
View at: Google Scholar
R. Torices and A. J. Muñoz-Pajares, “PHENIX: an R package to estimate a size-controlled phenotypic integration index,” Applications in Plant Sciences, vol. 3, no. 5, article 1400104, 2015.
View at: Publisher Site | Google Scholar
A. Dah, V. Hore, V. Iotchkova, and J. Marchini, “Network inference in matrix-variate Gaussian models with non-independent noise,” 2013, https://arxiv.org/abs/1312.1622.
View at: Google Scholar

Copyright

Copyright © 2021 Li-Chu Chien. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

254

Downloads

926

Citations