基于cis-SNP的多位点关联性分析方法比较及其在基因表达中的应用*

徐州医科大学公共卫生学院流行病与卫生统计学教研室(221004)

王鸣岐 安书仪 黄林玉 曾 平△ 黄水平△

【提 要】 目的 比较不同SNP集合方法检验常见变异和罕见变异与基因表达的关联,期望找到高效和稳健的统计方法。方法 对比Burden、SKAT、SKAT-O、MiST、ReLRT和GGRF等方法,通过GEUVADIS数据基因型和基因表达数据建立模拟研究,比较不同方法的一类错误率、统计效能及运算时间。基于常见和罕见cis-SNP,进一步通过GEUVADIS数据比较各种方法鉴别表达基因的能力。结果 模拟显示各种方法在不同变异位点个数下均能很好控制第一类错误，其中ReLRT和MiST的统计效能最高，MiST运算时间最短。在GEUVADIS数据分析中,ReLRT找到的表达基因最多。结论 各种检验方法基本都能有效控制第一类错误。模拟和实例数据显示ReLRT具有更高的检验效能,但计算速度较慢。

过去十年,全基因组关联研究(genome-wide association study,GWAS)发现了上万个与复杂疾病或性状相关联的单核苷酸多态性(single nucleotide polymorphism,SNP)位点[1]。然而,约90%的关联位点位于非编码区域,功能未知;理解关联位点的功能是目前遗传学研究面临的主要挑战之一[2]。GWAS和表达数量性状基因座(expression quantitative trait loci,eQTL)研究显示,许多关联SNP和eQTL是一致的[3]，即关联SNP不但影响复杂表型,还影响基因表达。这暗示关联SNP可能通过调控基因表达对复杂表型发挥作用。在eQTL研究中发现受SNP调控的基因至关重要,这种基因称为eGene[4];识别eGene为进一步研究SNP功能和理解疾遗传致病机理提供了重要的观点[4]。

单位点分析是常用的关联性检验方法,简单但统计效能低。相反,同时分析多个位点与复杂表型之间的关联可提高统计效能,称为SNP集合方法,在GWAS中获得了广泛引用[5]。然而,SNP集合关联方法在基因表达数据中的相对表现尚不清楚。与复杂疾病的关联研究一样,影响基因表达的SNP位点等位基因频率(minor allele frequency,MAF)分布也十分广泛,可能既有常见变异位点,也有罕见变异位点。本文中将MAF≥1%的SNP定义为常见变异,将MAF<1%的SNP定义为罕见变异[6]。我们的主要目的在于对比和评价SNP集合方法在eGene检验中的应用,在GEUVADIS数据[7]中通过分析常见变异和罕见变异与基因表达的关联,期望找到一个高效的、稳健的统计方法来检测eGene,也为进一步深入理解疾病遗传基础提供指导。

资料和方法

1.数据和方法

GEUVADIS数据包含五个不同人群 465个个体的基因表达量[7],基因型来源于千人基因组项目[8]。通过质量控制后保留15810个基因和7072917个SNP[9],按照MAF>1%将所有基因的cis-SNP分为常见变异和罕见变异,分别进行关联研究。本文主要比较以下方法:序列核关联检验(sequence kernel association test,SKAT[10])、优化序列核关联检验(sequence kernel association test-optimal,SKAT-O)[11]、负荷检验(burden test)[12]、混合效应得分检验(mixed-effects score test,MiST)[13]、限制性似然比检验(restricted likelihood ratio test,ReLRT)[14,15]和广义遗传随机场方法(generalized genetic random field method,GGRF)[16-17]。按照原理不同,上述方法大致可分为两类:将一组SNP位点合并为单个位点后进行关联分析(如Burden)和基于方差成分[18]的尺度参数检验(如SKAT和ReLRT[19])。由于每个遗传变异对其共同关联作用的贡献大小不同,需要考虑权重来平衡其中的差异性问题。对常见变异选择权重ωj=Beta(MAFj,1,1)[20],对罕见变异选择权重ωj=Beta(MAFj,1,25)[21];其中,MAFj指第j个SNP的最小等位基因频率,Beta指Beta分布的密度函数值。

SKAT通过得分函数(score function)检验随机效应的方差分量来研究罕见变异与疾病之间的关联,SKAT-O在SKAT检验中增加反映罕见变异关系结构的相关矩阵,该相关矩阵包含一个参数ρ,当ρ=1时SKAT-O退化为负荷检验,当ρ=0时SKAT-O退化为一般的SKAT;负荷检验的基本思想是首先将罕见变异在一定功能性区域内作为一个变量,然后对这些变量去进行检验;MiST方法,其思想是采用分层建模,利用已知的变异特征和利用位点信息来增强识别相关变异的能力;限制似然比检验是基于一组罕见变异可以被看作随机效应的线性混合效应模型而提出的,针对罕见变异与复杂疾病之间关联性分析的检验方法,需要对原模型与备择模型都要做出估计;GGRF方法认为,如果遗传变异是全部都与表型相关联的,个体间的基因相似性会导致他们的表型相似性。在GGRF方法中,如果把每个个体映射到欧几里得空间里,每个个体都有一个被其测序基因决定的位置。如果基因型与表型有关联性,则预计欧几里得空间中相邻的个体比距离远的个体有更多相似的表型。基于这一理念,GGRF模拟把每个个体的表型平均值作为剩余个体的表型加权平方和的一个线性函数,其权重取决于个体的基因相似性。

2.模拟设置

从GEUVADIS数据中抽取连续的遗传片段获得SNP的基因型G,通过Y=∑βjGj+ε产生基因表达,其中,βj是SNP效应,服从均数为0的正态分布,方差选择不同的数值;ε服从均数为0,方差为1的正态分布。设定基因型G包含10、30、50、80和100个SNP位点,检验水准为α=0.001。分别评价Burden、SKAT、SKAT-O、MiST、ReLRT和GGRF的第一类错误(模拟50000次)、统计效能(模拟1000次)和运算时间。本文分析均在R软件中完成,Burden、SKAT和SKAT-O通过软件包SKAT[21]执行,MiST方法通过软件包MiST[13]执行,ReLRT通过ReLRT函数[22]执行,GGRF通过GGRF函数[16-17]执行。

结 果

1.模拟数据结果

(1)第一类错误率

如图1A所示,在检验水准0.001时,在常见变异位点数为30、50、80以及100个时,犯第一类错误率最高的是MiST,其次分别是Burden、ReLRT、SKAT-O三种方法；变异位点数为10时,第一类错误率最高的是SKAT-O和ReLRT,在该检验水准下,对于常见变异的检验,犯第一类错误率最低的是GGRF和SKAT。从图1B可以看出,当罕见变异位点数为10、30时,犯第一错误率最高的是Burden和SKAT-O,当罕见变异位点数为50时,GGRE和ReLRT的第一错误率最高,当罕见变异位点数为80和100时,Burden和ReLRT的第一错误率最高,在该检验水准下,第一错误率最低的是SKAT和MiST。

(2)统计效能

如图1C所示,当常见变异位点数较少时(如10和30),统计效能最高的是ReLRT,其次是SKAT、SKAT-O和MiST,最低的是Burden,当变异位点增至50、80和100时,ReLRT和GGRF获得最高的统计效能,其次是SKAT、SKAT-O和MiST,最低的是Burden。由图1D可看出,MiST在各位点数情况下获得最高统计效能,而Burden获得最低的统计效能。

(3)运算时间

如表1所示,对于常见变异,MiST的运算时间在现有方法中最短,其次依次为Burden,SKAT,SKAT-O和ReLRT,GGRF的运算时间在现有方法中最长。而对罕见变异来说,MiST的运算时间在现有方法中最短,其次依次为Burden,SKAT,SKAT-O和GGRF,ReLRT的运算时间在现有方法中最长。

2.GEUVADIS 数据分析的结果

(1)通过常见变异检验到的eGene个数

在15,810个基因中,检验水准为0.05时,ReLRT找到的eGene最多,其次为SKAT;检验水准为0.01及0.001时,ReLRT找到的显著基因最多,分别为3681及2680个,其次为GGRF;Burden找到的显著基因最少,分别为1904、996及493个。结果见表2。

A:针对常见变异的不同统计方法的第一类错误率;B:针对罕见变异的不同统计方法的第一类错误率;C:针对常见变异的不同统计方法的统计效能;D:针对罕见变异的不同统计方法的统计效能。

图1 不同统计方法针对常见和罕见变异的在不同变异位点数下的第一类错误率和统计效能比较(检验水准α=0.001)

表2 通过常见变异检验到的eGene个数

(2)通过罕见变异检验到的eGene个数

通过罕见变异,在15,810个基因中检验水准为0.05时,GGRF找到的eGene最多,有2190个,其次为ReLRT;检验水准为0.01及0.001时,ReLRT找到的eGene最多,分别为1013及495个,其次为GGRF;Burden找到的eGene最少,分别为793、253及76个。结果见表3。

表3 通过罕见变异检验到的eGene个数

(3)常见变异与罕见变异表型之间的运算时间

不同方法对于GEUVADIS数据中常见变异与基因表达分析中,MiST的运算时间为0.5小时,在现有方法中运算时间最短,其次为Burden、SKAT、SKAT-O、ReLRT。GGRF的运算时间为18.4小时,在现有方法中运算时间最长;不同方法对于GEUVADIS数据中罕见变异与基因表达分析中,GGRF的运算时间为3.1小时,在现有方法中运算时间最短,其次为MiST、Burden、SKAT和SKAT-O。ReLRT的运算时间为598.8小时,在现有方法中运算时间最长。结果见表4。

表4 GEUVADIS数据常见变异与罕见变异与基因表达分析的运算时间(小时)

讨 论

本研究发现，Burden,SKAT,SKAT-O,Mist,ReLRT和GGRF六种方法在不同的检验水准、不同的变异位点数下,对第一类错误的控制都表现良好;但对于常见变异和罕见变异而言,表现各不相同。在常见变异中,第一类错误率最高的为GGRF和MiST,然而SKAT的第一类错误率在各种条件下均保持很低的水平。在罕见变异中,SKAT和MiST的表现仍然保守,第一类错误率处于低水平。MiST与SKAT相似,在本质上属于一种基于得分检验的关联分析方法,故具有一定的保守性,在一定样本量的保证下,基于方差成分检验的方法对于第一类错误率的控制基本不受检查位点个数影响,这也是随机效应检验优于固定效应检验之处。然而这种保守性会对检验效能产生负面的影响。

各方法的统计效能与变异位点数有关。在常见变异中,统计效能最高的一般都是ReLRT,最低的则是Burden,而在罕见变异中,MiST方法则表现出更好的统计效能,Burden方法最差。因此,针对常见变异和罕见变异位点的关联性分析,ReLRT和MiST方法各有优势。在运算速度上,可以看出,无论是对模拟数据的运算还是GEUVADIS实例数据进行运算,MiST均能表现出较快的运算速度,ReLRT在罕见变异与表型之间的检验的第一类错误率及统计效能的运算时间最长。对于常见变异,GGRF的运算时间最长,而在关联性分析中,ReLRT所需的时间最长,发现的显著基因最多。六种方法之间的应用情况不尽相同,其中负荷检验在变异的效应方向较为一致的时候应用可以获得最好的效果,然而这种情况并不常见;SKAT,SKAT-O和MiST方法在变异与性状关联性不明确或者仅有一小部分的关联性变异时,其检验比较稳健;GGRF与SKAT方法相比类似,主要区别在于所使用的检验方法与参数的不同,例如SKAT使用的是得分检验,而GGRF使用的是Wald类型的检验;ReLRT方法则比较全面稳健,但是其不足之处在于运算速度;因此,根据不同的数据和目的,需要我们对不同的方法进行灵活的选择。