多元表型与基因型的全基因组关联研究中的统计方法

勾建伟,刘应安,夏业茂



多元表型与基因型的全基因组关联研究中的统计方法

勾建伟,刘应安,夏业茂

南京林业大学理学院, 江苏 南京 210037

复杂疾病往往需要多元表型变量共同刻画，然而，常规的全基因组关联研究仅是检验单个表型，该策略往往由于忽略某些表型的信息而降低检验复杂疾病与基因关联性的效能。多元表型与基因型的全基因关联研究可以更有效揭示复杂疾病与基因的关联，同时对现有的统计方法提出挑战。本文综述了现有的处理多元表型与基因型的全基因组关联的统计方法，讨论了多元因变量多重回归模型框架下的回归系数与方差逆矩阵同时选择的组合稀疏方法的前景。该组合稀疏方法既能充分利用多元表型变量的相关性信息，增强变量选择的准确性，又能提供基因变异和表型关联程度的度量指标。同时基于多元因变量的稀疏方法可以推广到基因组学数据的整合分析中。

多元表型; 基因型; 全基因组关联性;多元统计分析

传统的全基因组关联研究(Genome-wide association studies，GWAS)是通过考察单核苷酸多态性(Single nucleotide polymorphism, SNP)位点与单个表型变量之间的关系。然而，刻画复杂疾病的量往往是错综复杂的，绝非仅仅单个表型变量就能够充分描述。因此，简单的用单一测量表型指标定义复杂疾病往往缺乏代表性。例如，身体质量指数(Body Mass Index，BMI)并不能全面地代表心脑血管疾病相关的肥胖特征，因为它只是粗略地测量了在给定体表面积的平均体重，并没有展示脂肪的分布。研究表明BMI、腰围、臀围这三个表型决定的体型更能反映心脑血管病风险[1]。因此，基于BMI、腰围、臀围的三个表型变量特征的GWAS与传统基于单个表型变量BMI的GWAS相比更合理。另一方面，从数量遗传学的观点看，同一种复杂疾病表型变量不仅受多个不同基因的影响，同一基因的变异也不只影响到一个表型变量，往往会同时影响到多个表型变量，这些表型变量之间往往具有相关性。可见，复杂疾病往往是由多个表型变量共同表现出来，同时检验多个表型变量比检验一个表型变量来分析疾病的相关性的效能要高。然而在传统GWAS背景下，通常将多元表型变量信息综合为一个得分(Score)单变量来进行基因关联分析。该策略往往会忽略某些对疾病发生、发展、转归有影响的表型信息，同时表型变量数据可能是二分类数据、有序数据、或者连续型数据以及三者间的任意混合，此时如何构造得分变量也面临挑战。传统的GWAS虽然已经取得丰硕成果，但距离完全揭开疾病的遗传密码相差甚远，已有研究[2]提出表型变量的复杂性、测量偏差及表型变量选取的不恰当都会大大削弱基因的关联性研究，造成信息损失。基于单个表型与基因的关联性分析难以全面揭示遗传与复杂疾病的联系。因此，我们需要通过多元表型与基因型的全基因关联研究，充分利用众多相关表型变量，探测出与疾病相关联的基因变异位点，从而为疾病的发病机制研究提供依据。

1 处理多元表型与基因型的全基因组关联研究的现有统计方法

传统的检验多元表型与基因相关性的方法是逐个对表型单变量与基因关联性进行检验，然后采用多重检验矫正。然而同一疾病的多元表型变量往往具有相关性，多重检验校正会降低效能。已有的研究利用多元表型变量间的相关性，从多元表型变量中提取类似潜变量因子的策略进行降维，然后再进行关联性分析，现有的统计方法概括如下。

1.1 多元因变量的主成分分析方法

通过对因变量进行主成分析，消除因变量之间的多重相关性[3]。对多表型的GWAS中，大多数研究[4,5]是采用前几个因变量的主成分与基因位点进行线性回归分析。然而，对于选取主成分的个数没有严格标准。已有研究[6]通过模拟不同相关性的表型变量，比较了不同情形下的多元因变量的主成分分析方法，并得出仅采用少数几个因变量主成分进行相关分析往往会降低检验的效能这一结论。同时，这类多元因变量的主成分方法仅仅综合了因变量的信息，对于超高维的基因自变量并没有给出有效的降维分析策略。

1.2 典型相关分析方法

分别提取自变量集与因变量集的最大主成分，通过两个主成分的相关关系推测自变量集与因变量集之间的相关关系[7]。由于典型相关分析也采用了主成分提取的思路，因此典型相关成分也拥有主成分降维的性质。已有研究[8-10]采用典型相关分析，对多表型变量与单个基因位点或者多个基因位点的相关性进行检验。然而，典型相关分析只是通过优化样本之间的协方差准则来求解典型变量，并没有考虑原始样本在高维空间中分布的结构，当GWAS这种相对于变量数目来说样本无法充分大的数据结构下，往往受样本影响，不太稳定。

1.3 偏最小二乘回归方法

偏最小二乘方法是多元因变量对多自变量的回归分析方法[11]。该方法考虑到自变量与因变量的相关程度，确保提取的成分对因变量解释能力强。已有研究[12]利用偏最小二乘回归方法对多元表型变量与多个位点集的相关性进行检验。然而，该方法也存在一些缺点，无法给出模型参数检验统计量的分布，只能采用置换检验(Permutation Test)得到检验的值。置换检验运算对于高维数据运算缓慢，无法应用在超高维的全基因组水平上，只能用于后期的验证性研究。同时，偏最小二乘提取出潜变量后，仅执行假设检验而不能参数估计，因此无法给出基因位点与表型变量相关性大小的度量，从而无法评价基因变异与疾病相关性大小。

1.4 非参数检验方法

若干用于多元表型GWAS的非参数检验方法也被提出，Zhang等[13]提出了广义Kendall’s秩检验用于检验单个基因位点与多表型变量的关联性。朱文圣等[14]提出改进的广义Kendall’s秩检验，研究了包含协变量调整下的基因与多表型变量的相关性检验。Van der Sluis等[15]提出组合值法，对每个表型变量关联分析的值提出一个组合值，可以有效的处理多元表型与基因关联性分析。

1.5 其它多元因变量的方法

基本是上述所有方法的不同扩展或改进，上面提到的方法主要是将基因位点作为一个集合(SNPs集)，检验SNPs集与元多表型变量的相关性，适用于变量维数和样本数相差不大的情形，然后通过多重检验应用到全基因组。对于高维自变量数据结构的全基因关联研究，另外一条思路是降维。通常的做法就是在上面的方法上增加惩罚函数，提出对应的稀疏方法。例如，已有研究[16]提出偏最小二乘方法与1惩罚结合的稀疏偏最小二乘方法，对多元表型变量与基因位点集的相关性给出了检验。该方法被应用到多元表型的若干GWAS中[17,18]。

上述所有方法不预先设定模型结构形式，所选方法都由所观测到的样本数据所决定，因此具有较大的适应性，在降低建模偏差方面有较大优势，但也存在明显的局限和不足：首先，变量的维数较高时，上述方法所涉及的方差及协方差往往是奇异的，给数值计算与理论性质带来了难度。其次，上述方法对基因变异的相关性都是采用的检验，给出检验的值，但很难给出基因关联效应的估计值，无法合理解释基因变异影响性状表型的大小。从理论和实际应用角度讲，好的多元表型与基因关联分析方法应该具备以下几个特点：能够给出更快速更有效的估计值，提供度量基因变异和表型变量关联大小的客观指标；能够有效给出表型变量与基因位点回归参数的渐进分布；能够对遗传数据的潜在模型提出更好的解释，揭示疾病和基因变异之间的复杂关联性。

2 稀疏多元因变量多重回归方法的前景和意义

稀疏回归或者惩罚回归模型是在损失函数的基础上增加惩罚函数，通过调整参数的选择，将不显著参数压缩为零，对显著参数进行很小压缩或者不压缩，并且给出该参数的估计。稀疏回归具有以下优势：允许变量的维数超过样本总数，通过变量选择选出显著的变量；能够在进行变量选择的同时，给出模型参数的优良估计。从而，稀疏回归方法在高维数据分析中有着良好的表现，是近年来高维数据领域中十分热门的课题，比如：LASSO[19]、LARS[20]、Elastic net[21]、SCAD[22,23]、自适应LASSO[24]、成组LASSO[25]、自适应成组LASSO[26]、两水平惩罚方法[27]以及图LASSO方法[28]。

另一方面，多元因变量多重线性回归模型建立了多元因变量与多自变量之间如下的线性关系：=+其中表示×的因变量矩阵，表示×的设计矩阵，是×的未知系数矩阵，表示×的误差量矩阵，并假定的个行向量独立同分布于元正态分布(0;Ʃ)。

其中Ω=Σ-1，称为精度矩阵，刻画了元因变量的两两条件相关性。

模型中待估的回归参数刻画预测变量与因变量的线性关系，精度矩阵表示多元因变量的两两条件相关性。当变量维数特别大时，假设仅有一部分自变量与因变量相关，其余变量都是噪声变量，此时系数矩阵应该是稀疏的；当因变量维数特别大，若干因变量之间的相关性应该比较小，此时精度矩阵Ω应该是稀疏的。当回归系数矩阵和精度矩阵都满足稀疏性时，可以通过参数的惩罚函数，实现参数的稀疏性约束。稀疏的多元因变量多重线性回归的模型估计转变为：

同时，稀疏多元因变量多重回归模型不仅可以实现变量选择，还可以给出回归系数矩阵和精度矩阵的估计。回归系数的估计值可以反映自变量和因变量的相关性；利用精度矩阵与高斯图解模型的关系，我们可以根据精度矩阵刻画出因变量的相关关系图结构，进一步深入理解复杂疾病的遗传机制。

综上，在多元因变量多重线性回归框架下发展和应用惩罚回归方法，将同时具备上面提到的三个特点。稀疏的多元因变量多重线性回归将为基于多元表型与基因型的复杂疾病关联研究提供更加灵敏和强有力的统计分析工具。

3 讨论

针对多元表型与基因型数据的复杂疾病关联分析，借助组合稀疏回归模型，提出在多元因变量多重回归的框架下进行组合惩罚回归的设想，根据不同的数据结构选择和构造不同的惩罚函数1与2。稀疏的多元因变量多重回归模型中回归系数矩阵和精度矩阵Ω同时约束的最优解求解往往不易求。现有的方法主要是分为两类：令2(Ω)=0的降秩方法与令1()=0的高维协方差选择的图LASSO(graph LASSO)方法。这两类方法无法同时考虑自变量和因变量进行变量选择的高维数据结构。我们可以借助惩罚似然与惩罚条件似然的关系，探索同时约束回归系数矩阵和精度矩阵Ω的最优解算法。

针对现有多元表型关联性分析检验方法统计量分布未知和实际应用受限等问题，需要给出稀疏回归模型参数估计理论框架，包括研究估计值的无偏性、相合性以及推导渐进分布；结合坐标下降优化算法，给出参数计算快速有效的算法。方法比较、估计值的理论研究（无偏性，相合性以及渐进分布）、数值模拟和实例分析等方面系统探索新提出的方法，为多元表型的GWAS数据和基因组学数据的整合分析提供实用、高效和灵敏的关联性分析工具其顺利实施不但能够丰富和发展关联性分析理论，也对进一步深入理解疾病遗传基础和解释遗传缺失具有重要意义。

[1] Wells JC, Treleaven P, Cole TJ. BMI compared with 3-dimensional body shape: the UK National Sizing Survey[J]. Am J. Clin. Nutr., 2007(85):419-425

[2] Sluis SVD, Verhage M, Posthuma D,. Phenotypic complexity, measurement bias, and poor phenotypic resolution contribute to the missing heritability problem in genetic association studies[J]. Plos One, 2010,5(11):e13929

[3] Cadima J, Jolliffe IT. Loading and correlations in the interpretation of principle components[J]. Journal of Applied Statistics, 1995(22):203-214

[4] Liu F, Van DLF, Schurmann C,. A genome-wide association study identifies five loci influencing facial morphology in Europeans[J]. Plos genetics, 2012,8(9):e1002932

[5] Zhang F, Guo X, Wu S,. Genome-wide pathway association studies of multiple correlated quantitative phenotypes using principle component analyses[J]. Plos One, 2012,7(12):e53320

[6] Aschard H, Vilhjálmsson BJ, Greliche N,. Maximizing the power of principal-component analysis of correlated phenotypes in genome-wide association studies[J]. The American Journal of Human Genetics, 2014,94(5):662-676

[7] Hotelling H. Relations between two sets of variates[J]. Biometrika, 1936,28(3/4):321-377

[8] Galesloot TE, Van SK, Kiemeney LA,. A comparison of multivariate genome-wide association methods[J]. Plos One, 2014,9(4):e95923

[9] Seoane JA, Campbell C, Day INM,. Canonical Correlation Analysis for Gene-Based Pleiotropy Discovery[J]. Plos Computational Biology, 2014,10(10):e1003876

[10] Tang CS, Ferreira MAR. A gene-based test of association using canonical correlation analysis[J]. Bioinformatics, 2012,28(6):845-850

[11] Geladi P, Kowalski BR. Partial least-squares regression: a tutorial[J]. Analytica chimica acta, 1985,185(86):1-17

[12] Xue F, Li S, Luan J,. A latent variable partial least squares path modeling approach to regional association and polygenic effect with applications to a human obesity study[J]. Plos One, 2012,7(2):e31927

[13] Zhang H, Liu CT, Wang X. An association test for multiple traits based on the generalized Kendall’s tau[J]. Journal of the American Statistical Association, 2010,105(490):473-481

[14] Zhu W, Jiang Y, Zhang H. Nonparametric covariate-adjusted association tests based on the generalized Kendall's Tau.[J]. Journal of the American Statistical Association, 2012,107(497):1-11

[15] Van dSS, Posthuma D, Dolan CV. TATES: efficient multivariate genotype-phenotype analysis for genome-wide association studies[J]. Plos genetics, 2013,9(1):e1003235

[16] Chun H, Ballard DH, Cho J,Identification of association between disease and multiple markers via sparse partial least‐squares regression[J]. Genetic epidemiology, 2011,35(6):479-486

[17] Allen GI, Peterson C, Vannucci M,. Regularized partial least squares with an application to NMR spectroscopy[J]. Statistical Analysis and Data Mining the ASA Data Science Journal, 2013,6(4):302-314

[18] Rinnan A, Andersson M, Ridder C,. Recursive weighted partial least squares (rPLS): an efficient variable selection method using PLS[J]. Journal of Chemometrics, 2014,28(5):439-447

[19] Tibshirani R. Regression shrinkage and selection via the lasso[J]. Journal of the Royal Statistical Society: Series B (Methodological), 1996,58(1):267-288

[20] Efron B, Hastie T, Johnstone I, et al. Least angle regression[J]. Institute of Mathematical Statistics, 2004,32(2):407-451

[21] Zou H, Hastie T. Regularization and variable selection via the elastic net[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2005,67(2):301-320

[22] Fan J, Li R. Variable Selection via Nonconcave Penalized Likelihood and its Oracle Properties[J]. Journal of the American Statistical Association, 2001,96(456):1348-1360

[23] Kim Y, Choi H, Oh HS. Smoothly clipped absolute deviation on high dimensions[J]. Journal of the American Statistical Association, 2008,103(484):1665-1673

[24] Zou H. The Adaptive Lasso and Its Oracle Properties[J]. Journal of the American Statistical Association, 2006,101(476):1418-1429

[25] Yuan M, Lin Y. Model selection and estimation in regression with grouped variables[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2006,68(1):49-67

[26] Wang H, Leng C. A note on adaptive group lasso[J]. Computational Statistics & Data Analysis, 2008,52(12):5277-5286

[27] Seetharaman I. Consistent bi-level variable selection via composite group bridge penalized regression[D].Kansas, USA: Kansas State Univesity, 2013

[28] Fang Y, Wang R, Dai B,. Graph-based learning via auto-grouped sparse regularization and kernelized extension[J]. Transactions on Knowledge & Data Engineering, 2015,27(1):142-154

[29] Yuan M. Dimension reduction and coefficient estimation in multivariate linear regression[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2007,69(3):329-346

Statistical Method of Multiple Phenotype and Genotype in Genome-wide Association Study

GOU Jian-wei, LIU Ying-an, XIA Ye-mao

210037,

The complex disease is often identified with multiple traits. However, the common genome-wide association studies only test the association between a single phenotype trait and a large number of single nucleotide polymorphisms (SNPs), Which will lead to loss in statistical power when multivariate phenotype are involved. The genome-wide association study of multiple phenotype and genotype can more effectively reveal the association between complex diseases and genes, and challenge existing statistical methods. This paper reviewed the existing statistical methods for dealing with genome-wide associations between multivariate phenotype and genotype, and discussed the prospects of combined sparse methods for simultaneous selection of regression coefficients and variance inverse matrices under the framework of multivariate dependent multivariate regression models. The combined sparse method not only improved the accuracy of variable selection by utilizing joint information among multivariate phenotype, but also offered a measure of the association between genotype and phenotype. At the same time, the method could extend the proposed sparse multivariate regression methods to the integrative analysis of genomic data.

Multiple phenotype; genotype; genome-wide association; multivariate statistical analysis

O212

A

1000-2324(2018)05-0906-05

10.3969/j.issn.1000-2324.2018.05.038

2017-09-02

2017-10-13

国家自然科学基金(11471161);江苏省高校基金(15KJB110010)

勾建伟(1982-),女,博士,讲师,研究方向为统计基因组学. E-mail:gjw1983@139.com