利用全基因组表达数据对比评价九种遗传预测模型方法*

徐州医科大学公共卫生学院流行病与卫生统计学系(221004) 余星皓 刘 兵 曾 平 黄水平

【提 要】 目的 研究稀疏模型(Lasso、ENET、ssLasso、贝叶斯变量选择回归模型(BVSR))与多基因模型(线性混合模型(LMM)、贝叶斯稀疏线性混合模型(BSLMM)、狄利克雷回归模型(DPR))等九种遗传预测方法在全基因组表达数据中对复杂疾病的遗传预测表现。方法 通过模拟研究评价每种方法在不同的较大基因稀疏程度和不同的遗传度下的预测精度,利用乳腺癌数据进行表型预测。结果 模拟结果显示预测方法在满足各自的模型假设时表现结果最好。在相同模拟假设情况下,随着遗传度的增高,模型的预测准确性也逐渐增高。BVSR运算速度和BSLMM运算速度相似,由于迭代次数的影响,BVSR与BSLMM的运算速度低于LMM。实际的乳腺癌数据显示BSLMM和DPR的预测精度优于其他方法。结论 BSLMM和DPR在不同模拟情形下和真实数据中均表现出稳健的预测能力,值得在实际应用中推荐。

全基因组关联研究(genome-wide association study,GWAS)已经发现了上万个单核苷酸多态性(single-nucleotide polymorphism,SNP)位点与数百种复杂疾病存在关联,为表型变异的遗传基础提供了前所未有的视角[1-7]。GWAS的成功也极大促进了利用遗传信息(除环境和生活方式信息外)进行复杂疾病风险预测评估的研究[8]。在动物或植物中,具有遗传标记的准确表型预测可以帮助选择理想的育种个体,从而提高育种计划的有效性[9];在人群中,利用遗传标记的准确表型预测有利于复杂疾病的早期预防和干预、促进开发个性化药物,从而制定治疗方案和预测疗效[10]。此外,在跨组学数据分析中,表型预测也被认为是将功能基因组测序研究与GWAS相结合的关键步骤,可以在GWAS中构建更高效和具有可解释性的基因集测试,通过设置SNP权重,将其从预测模型中推断出来,并用于表达数量性状位点映射研究[11]。

研究者们发明了许多新的表型预测模型方法以利用全基因组遗传位点信息提高预测精度。表型预测模型本质上是对基因效应大小分布的假设不同,例如Henderson等人[12]提出的线性混合模型(liner mixed model,LMM),Speed等人[13]提出的最佳无偏估计模型(multiple best linear unbiased prediction,MultiBLUP),Tang等人[14]提出的Bayes Lasso模型。上述方法多属于参数模型,所假设的遗传效应分布是高度参数化的、由少数几个参数决定,其局限在于仅对特定满足模型假设的疾病预测效果良好,对其它不满足假设的疾病预测精度不佳。理论上,模型预测精度取决于遗传效应先验与真实分布的接近程度。然而,真实的遗传结构往往未知,由于疾病的不同，在遗传度、关联位点个数、最小等位基因频率和效应大小等方面存在诸多差别。另外,这些方法先前主要利用全基因组SNP位点信息进行复杂疾病的预测;而越来越多的研究表明基因表达在复杂疾病的发生、发展和转归中起着十分重要的作用,是SNP对复杂疾病起作用的重要分子中介。因此,通过全基因组表达数据预测复杂疾病具有重要的科学意义。通过基因表达预测复杂疾病与利用SNP预测复杂疾病的差别在于，基因表达数据的变量较少,效应可能更强,作为一个整体的功能单位,基因相对SNP更具有结果的可解释性。那些通过全基因组SNP位点获得复杂疾病预测结论未必适合基因表达预测。

本文主要对九种遗传预测模型方法进行比较,这些方法大致可分为两类:稀疏模型和多基因模型。前者包括Lasso、elastic net(ENET)、spike-and-slab Lasso GLMs(ssLasso)、贝叶斯变量选择回归模型(Bayesian variable selection regression,BVSR);后者包括线性混合模型(LMM)、贝叶斯稀疏线性混合模型(Bayesian sparse linear mixed model,BSLMM)、狄利克雷回归模型(Dirichlet process regression,DPR)等。本文通过癌症患者基因表达数据模拟各种复杂疾病可能的遗传结构,比较上述方法在不同情况下对于复杂疾病的预测能力，进一步通过比较各方法在真实肿瘤数据中的预测能力评价这些方法在基因表达数据中的表型预测精度。

材料与方法

1.方法

比较九种模型(稀疏模型:BVSR[15]、Lasso[16]、ENET[17]和ssLasso,多基因模型:LMM、Linear-BSLMM、Probit-BSLMM[18]、DPR.VB和DPR.MCMC[19])在遗传预测中的能力,这些模型方法被广泛运用于人类复杂疾病的表型预测以及动植物繁殖的基因选择中。具体而言,BVSR假设较少部分的基因变异对复杂疾病或性状产生效应,单个基因变异对表型的效应较大,假设总体效应大小符合正态分布和零点分布的混合分布;Lasso是在模型的离均差平方和基础上增加一个绝对值的惩罚项,使其总和小于等于一个常数,是一种压缩估计方法;ENET是一种Lasso与岭回归组合后的回归分析,在绝对值惩罚项基础上同时施加平方和惩罚项;ssLasso假设基因的系数服从混合双指数先验,根据不同的基因型与表型数据,产生适合的收缩系数,并保留系数较大的基因;LMM假设所有的基因对表型均产生效应,且效应均较小,同时服从正态分布;BSLMM假设基因对表型的效应大小服从两个正态分布组成的混合分布,两个正态分布的混合程度决定于不同的混合比例,需要注意的是,针对病例对照研究(假设表型编码为0-1),BSLMM可直接将0-1当作连续变量通过线性模型处理,称为Linear-BSLMM,或将Probit链接函数通过广义线性模型处理,称为Probit-BSLMM,图1中分别以BSLMM1和BSLMM2表示;DPR假设基因效应大小服从一个无穷混合正态分布,属于非参数贝叶斯模型范畴,DPR的参数估计可通过传统的马尔科夫蒙特卡罗估计,称为DPR.MCMC,也可通过变分贝叶斯方法估计,称为DPR.VB。

LMM、BVSR、Linear-BSLMM和Probit-BSLMM方法采用GEMMA[18]软件执行,DPR.VB、DPR.MCMC等方法采用DPR[19]软件执行,Lasso、ENET通过glmnet软件包执行,ssLasso通过BhGLM[20]软件包执行。

2.模拟数据

为更加接近真实数据,使用916例乳腺癌患者的基因表达数据,模拟三种情境下的表型数据:(i)所有基因均对表型产生效应,分别服从正态分布、t分布或Laplace分布,设置遗传度分别为0.2、0.5或0.8;(ii)仅有小部分基因(分别为10、50或100)对表型产生效应,服从正态分布,大部分剩余的基因对表型均不产生效应,遗传度分别为0.2、0.5或0.8;(iii)所有基因对表型均产生效应,效应基因分为两部分:小部分效应较大的基因和较大部分效应较小的基因,两组基因的效应大小都服从正态分布,两部分占遗传度的比例分别为0.2和0.8;较大效应的基因数量分别为10、50或100;遗传度分别为0.2、0.5或0.8。每种情景分别重复模拟20次,每一次重复中随机将数据中90%的个体作为训练集,剩余的10%作为测试集。通过训练集数据估计每种方法的模型参数,并在测试集数据中比较不同方法的预测能力。通过表型与预测值的曲线下面积(area under the curve,AUC)评价预测效果。

3.真实数据

真实数据全部来源于癌症和肿瘤基因图谱(The Cancer Genome Atlas,TCGA),通过加利福尼亚大学基因组浏览器(UCSC Xena,https://xenabrowser.net/)下载乳腺癌(Breast Cancer,BRCA)数据。原始数据包括1215例患者临床数据和1219例患者的20530个基因表达数据,首先对两份数据进行合并,删除重复的患者和男性患者,同时删除2840个零表达值超过50%的基因,最终获得916例患者的17690个基因表达数据。在数据分析前,首先将癌症的基因表达数据分成两部分:90%的患者为训练集数据,剩余10%的患者为测试集数据,对数据集执行20次重复。通过AUC评价不同模型的预测能力。

结果

1.模拟数据结果

情景i中(图1-i)当遗传度为0.2时,效应大小的分布服从正态分布的情况下,各模型的预测能力接近。效应大小的分布服从t分布和Laplace分布的情况下,多基因模型的预测能力均高于稀疏模型,多基因模型的预测能力大小接近。当遗传度为0.5或0.8时,三种情况下多基因模型中五种模型的预测准确性相对接近,优于稀疏模型,BVSR预测能力最低。当遗传度为0.8时,各统计学模型的预测结果AUC均在0.7左右,高于遗传度为0.2和0.5时的预测结果。

情景ii中(图1-ii)当遗传度为0.2或0.5时,稀疏模型的预测能力优于多基因模型,但是,随着效应基因数量和遗传度增加,稀疏模型与多基因模型之间的差异逐渐缩小;当遗传度为0.8:效应基因数量为10的时候,BSLMM和DPR的预测能力较强,其余多基因模型精度稍低。随着效应基因数量增加,稀疏模型与多基因模型之间的差异越小,当效应基因数量为100时,稀疏模型的预测能力要低于多基因模型。

情景iii中(图1-iii)当遗传度为0.2时,三种情况下多基因模型的预测效果均高于稀疏模型,其中BVSR的预测能力最差,linear-BSLMM、probit-BSLMM和DPR.MCMC的预测能力较强且较稳定。当遗传度为0.5时,多基因模型的预测效果优于稀疏模型,BSLMM和DPR的预测能力较稳定,且精度较高。当较大效应基因数量增加时,Lasso和ENET的预测能力上升;当遗传度为0.8时,较大效应基因的数量不同时,模型之间预测结果相似。多基因模型的预测结果相似。稀疏模型的预测结果相对不稳。在不同情况下,多基因模型的预测能力均稍高于稀疏模型。各统计学模型的预测结果AUC均在0.75左右,高于遗传度为0.2或0.5时的预测结果。

图1 遗传度0.2(A)、0.5(B)、0.8(C)时三种模拟情景下Lasso(红)、ENET(绿)、ssLasso(深蓝)、BVSR(绿)、LMM(蓝)、BSLMM1(Linear-BSLMM)(黄)、BSLMM2(Probit-BSLMM)(紫)、DPR.VB(浅蓝)、DPR.MCMC(粉)之间AUC估计值的比较

2.真实数据结果

真实数据来源于上述的乳腺癌数据,将术后淋巴结阳性转移作为预测的表型数据(编码为0-1,视为连续变量应用线性模型处理)。经过计算,数据的真实遗传度为0.15。各方法表型预测结果AUC均在0.61左右(图2A),多基因模型的预测能力均高于稀疏模型。多基因模型中,五种方法的预测能力几乎无差异(虽然DPR和BSLMM轻微略高)。稀疏模型中,ENET的预测能力最强,Lasso与BVSR的预测能力接近,ssLasso的预测能力最低。同时,采用Logistic回归分析各基因的效应大小,并绘制相应的分布密度函数,与正态分布、Laplace分布和自由度为4的t分布相比较(图2B)。

3.运算时间

最后对比各种方法的运算速度,以乳腺癌数据为例。这里需要说明的是,贝叶斯预测方法的运算时间与所使用的迭代次数有关系,为了节省时间,这里只进行了100次迭代。结果显示(表1)，Lasso、ENET和ssLasso的运算速度最快,其次为LMM和DPR.VB;BVSR、BSLMM和DPR.MCMC运算速度相近,且时间均较长。

讨论

本文所提及的多种统计学方法在遗传统计学的多个方面都有较好的应用价值，例如估计遗传度[18]和关联映射[21]。本文主要致力于多种统计学方法在表型预测方面的研究。国内外已发表的文献[13,18-19,22-24]主要研究表型预测方法在GWAS中的比较,已取得较大成果,本文通过模拟不同的疾病遗传结构,比较了两类统计学方法在基因表达数据中的预测能力。模拟数据分析结果显示,各模型在表型预测上的表现不同。当模拟假设全部基因对表型均产生效应时,多基因模型的预测能力要高于稀疏模型;当模拟假设只有部分基因对表型产生效应时,稀疏模型的预测能力要高于多基因模型。当模拟假设小部分基因对表型产生较大的效应,其余大部分对表型产生较小的效应时,多基因模型中的BSLMM和DPR预测能力要高于其余模型。在相同模拟假设情况下,随着遗传度的增高,模型的预测准确性也逐渐增高。当假设全部基因对表型均产生效应、效应大小服从单个正态分布、遗传度值较小时,多基因模型与稀疏模型的预测能力相似。当遗传度较大时,多基因模型的准确性稍高于稀疏模型。当假设部分基因对表型产生效应,且效应大小服从单个正态分布时,随着相应基因数量的增加,稀疏模型的预测能力逐渐由高于多基因模型逐渐变化为与多基因模型预测能力相似,甚至低于多基因模型。真实数据分析结果显示不同统计学方法在表型预测的表现不同,各方法表型预测结果表型预测结果AUC均在0.61左右,多基因模型的预测能力均高于稀疏模型,其中BSLMM和DPR最优。由模拟研究和真实数据分析结果可以看出,模型预测的准确性与复杂疾病或性状的真实遗传结构有关,当复杂疾病的真实遗传结构与模型假设相吻合时,模型的预测表现较好,而当复杂疾病的真实遗传结构与模型假设差异较大时,模型的预测表现较差。这与以往研究的结论相同[13,18]。因此我们认为多基因的模型假设与复杂疾病的真实遗传结构相符合,少数较大效应的基因和大多数小效应基因共同对复杂疾病产生影响。

图2 BRCA九种方法预测能力的比较

表1 真实数据集中不同方法的运算时间(小时)

*:Lasso、ENET、ssLasso的调整参数均为通过100倍交叉验证选择的最佳值,BVSR、BSLMM、DPR.MCMC均为进行了100次MCMC迭代计算,括号中的值为标准差,均值和标准差都是基于模拟的20次重复计算。

本文研究的基因表达数据相对较小,三种方法的运算速度均较高,但是这些方法往往并不适用于较大数据集[22]。需要注意的是,Lasso、ENET和ssLasso的运算速度同时受到交叉验证次数的影响,本研究中只通过100倍交叉验证,并未进一步研究交叉验证数量对于运算速度的影响。LMM、BVSR和BSLMM使用GEMMA软件进行分析,同时适用于大规模数据集。BVSR和BSLMM的运算速度相似,由于迭代次数的影响,BVSR与BSLMM的运算速度低于LMM。有研究报道[18],当纳入模型的自变量较大时,BSLMM的运算速度将远高于BVSR。DPR通过DPR软件进行分析,其中DPR.VB与LMM运算时间接近,而DPR.MCMC运算时间同样受到迭代次数的影响,远高于其它方法。此外,计算时间还受到工作环境、计算机语言、以及基因和表型之间的关系等因素的影响。虽然目前对于表型预测方法的研究层出不穷,但是尚未有一种方法对于不同遗传结构的复杂疾病或性状都具有较高的预测能力且运算时间较快。同时,遗传度较低时,所有方法的预测表现均较差。此外,本文所纳入的方法均为通过基因表达数据进行复杂疾病或性状的预测,都未能考虑到在遗传数据中非常重要的分组结构。例如:SNP可以通过基因或者功能注释的形式分成不同的组;基因之间存在交互作用和共同作用,基因可以通过通路信息进行分组,功能相同的基因可以被划分到一组[25]。这些被忽视的遗传信息(分组结构、非线性效应和交互作用)可能会对预测结果造成一定的影响。已经有研究表明,充分利用多组学的信息可以提高遗传预测的精度[26],这也将是接下来的研究重点之一。