曾 平，趙 楊，陳 峰，金英良，張耀東

1)徐州醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室 江蘇徐州 221004 2)南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)系 南京 210009

病例-對(duì)照設(shè)計(jì)的似然比罕見變異關(guān)聯(lián)性檢驗(yàn)的構(gòu)建和模擬評(píng)價(jià)*

曾 平1)△，趙 楊2)，陳 峰2)，金英良1)，張耀東1)

1)徐州醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室 江蘇徐州 221004 2)南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)系 南京 210009

△男，1982年7月生，博士，講師，研究方向：高維數(shù)據(jù)分析和統(tǒng)計(jì)遺傳學(xué)方法研究，E-mail：zpstat@xzhmu.edu.cn

罕見變異；關(guān)聯(lián)性分析；logistic混合模型；似然比檢驗(yàn)

目的：發(fā)展適用于病例-對(duì)照設(shè)計(jì)的似然比罕見變異關(guān)聯(lián)性分析方法。方法：在logistic混合模型的框架下基于PQL算法建立偽數(shù)據(jù)，將二分類表型轉(zhuǎn)化為連續(xù)型表型的關(guān)聯(lián)性分析，然后借助線性混合模型的方差成分檢驗(yàn)執(zhí)行關(guān)聯(lián)性分析。采用Monte Carlo模擬評(píng)價(jià)該方法的有效性，并與現(xiàn)有方法進(jìn)行對(duì)比。結(jié)果：模擬顯示，在不同情況下包括似然比檢驗(yàn)在內(nèi)的所有統(tǒng)計(jì)檢驗(yàn)都能有效控制Ⅰ型錯(cuò)誤；在效應(yīng)方向相同情況下，Burden檢驗(yàn)、SKAT-O和MiST的統(tǒng)計(jì)效能高；在效應(yīng)方向不同的情況下似然比檢驗(yàn)優(yōu)于其他方法。結(jié)論：基于logistic混合模型和PQL算法的似然比檢驗(yàn)可有效用于病例-對(duì)照設(shè)計(jì)的罕見變異關(guān)聯(lián)性分析。

過去十年，盡管全基因組關(guān)聯(lián)研究取得了巨大成功，但對(duì)許多復(fù)雜疾病而言，已發(fā)現(xiàn)的絕大部分常見位點(diǎn)效應(yīng)都很弱，僅能解釋極少一部分遺傳度，未能取得預(yù)期的突破性進(jìn)展[1-3]。越來越多的證據(jù)表明罕見變異對(duì)復(fù)雜疾病有實(shí)質(zhì)性的貢獻(xiàn)[4]，支持的論據(jù)包括[5]：①進(jìn)化理論暗示與疾病相關(guān)的遺傳變異應(yīng)該是近期發(fā)生的并且是罕見的。②人群遺傳研究的許多經(jīng)驗(yàn)證據(jù)表明對(duì)人類健康有害的變異是罕見的。③一些心理疾病由罕見的拷貝數(shù)變異引起。④許多不常見的家族性疾病由罕見遺傳位點(diǎn)變異所致。⑤由罕見變異引起的虛擬關(guān)聯(lián)可解釋部分常見變異與疾病的關(guān)聯(lián)[6]。然而，由于最小等位基因頻率極低，很難對(duì)罕見變異進(jìn)行關(guān)聯(lián)性研究。作者主要針對(duì)病例-對(duì)照設(shè)計(jì)，在廣義混合效應(yīng)模型框架下構(gòu)建基于似然比檢驗(yàn)的罕見變異關(guān)聯(lián)性分析方法，通過Monte Carlo評(píng)價(jià)該方法的性質(zhì)、表現(xiàn)。

1 方法原理和模擬設(shè)置

1.1 方法原理 設(shè)Y=1和0分別表示病例和對(duì)照，X為協(xié)變量，G為K維基因型，樣本量為n。采用logistic混合效應(yīng)模型描述Y、X和G之間的關(guān)系：

式中，W=Beta(MAFk|1,25)，為K維的對(duì)角權(quán)重矩陣，Beta(a1,a2)為Beta函數(shù)，MAF表示最小等位基因頻率，τ2為方差參數(shù)。檢驗(yàn)G是否與Y存在關(guān)聯(lián)等價(jià)于檢驗(yàn)H0：τ2=0。采用懲罰擬似然算法(penalized quasi-likelihood，PQL)估計(jì)模型的參數(shù)，通過R軟件中的函數(shù)glmmPQL[7]。具體而言，PQL算法在偽對(duì)數(shù)極大(ml)或限制性極大(reml)似然函數(shù)和工作變量Y′間迭代：

lg|X′D-1X|]

Y′=Xβ+Gb+ε

其中R是對(duì)角元素為μ(1-μ)的n維對(duì)角矩陣。當(dāng)PQL算法收斂時(shí)，工作應(yīng)變量為

則

V=τ2G*WWG*′+In

這樣就將病例-對(duì)照設(shè)計(jì)的方差成分檢驗(yàn)轉(zhuǎn)化成為了連續(xù)型表型的方差成分檢驗(yàn)。建立(偽)似然比統(tǒng)計(jì)量

也即是，在執(zhí)行方差成分檢驗(yàn)時(shí)只要將偽數(shù)據(jù)當(dāng)作真實(shí)數(shù)據(jù)即可。

1.2模擬設(shè)置 采用COSI軟件模擬聯(lián)合模型產(chǎn)生歐洲人群的罕見遺傳變異基因型[8]。協(xié)變量分別為相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)變量(X1)和率為0.50的二項(xiàng)分布隨機(jī)變量(X2)。Y從二項(xiàng)分布產(chǎn)生

其中，C為關(guān)聯(lián)位點(diǎn)數(shù)，b為效應(yīng)；Ⅰ型錯(cuò)誤控制時(shí)b=0，統(tǒng)計(jì)效能評(píng)價(jià)時(shí)b=0.4|lgMAF|。Ⅰ型錯(cuò)誤控制和統(tǒng)計(jì)效能模擬次數(shù)分別為10 000和5 000，檢驗(yàn)水準(zhǔn)α設(shè)定為0.01，Ⅰ型錯(cuò)誤和統(tǒng)計(jì)效能由P≤α(=0.01)的比例估計(jì)。統(tǒng)計(jì)分析在R軟件中完成，Burden檢驗(yàn)、SKAT和SKAT-O執(zhí)行軟件包SKAT[9]，MiST[10]和GenRF[11]均執(zhí)行相應(yīng)函數(shù)。分別考察固定基因片段和不同基因片段條件下的Ⅰ型錯(cuò)誤控制和統(tǒng)計(jì)效能分析。

2 結(jié)果

2.1 模擬1：固定基因片段條件下的Ⅰ型錯(cuò)誤控制和統(tǒng)計(jì)效能分析 選擇30%的基因片段，樣本量為400、600、800和1 000時(shí)Ⅰ型錯(cuò)誤控制情況見表1。結(jié)果表明在固定基因片段和不同樣本量條件下，SKAT和MiST展現(xiàn)出保守的傾向，而其他方法基本上都能較好地控制Ⅰ型錯(cuò)誤。

表1 固定基因片段和不同樣本量條件下的Ⅰ型錯(cuò)誤控制

α=0.01。

選擇30%的基因片段(罕見變異平均個(gè)數(shù)約70)，將其中30%的罕見變異設(shè)為關(guān)聯(lián)位點(diǎn)，并假設(shè)3種情況：①所有關(guān)聯(lián)位點(diǎn)效應(yīng)方向?yàn)檎?；?0%關(guān)聯(lián)位點(diǎn)效應(yīng)方向?yàn)樨?fù)，70%效應(yīng)方向?yàn)檎?；?0%關(guān)聯(lián)位點(diǎn)效應(yīng)方向?yàn)樨?fù)，50%效應(yīng)方向?yàn)檎?。樣本量?00、600、800和1 000。為了解大樣本情況下似然比檢驗(yàn)的表現(xiàn)，可將樣本量增加到2 000。結(jié)果顯示，當(dāng)所有關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為正時(shí)，Burden檢驗(yàn)、SKAT-O和MiST的統(tǒng)計(jì)效能更高，其次是LRT和ReLRT，而SKAT、GenRF的效能最小(表2)。當(dāng)關(guān)聯(lián)位點(diǎn)效應(yīng)方向不一致時(shí)，LRT、ReLRT和SKAT表現(xiàn)出更高的效能，而GenRF、SKAT-O和MiST次之，Burden檢驗(yàn)的效能最低(表3、4)。相對(duì)效應(yīng)方向一致的情況，效應(yīng)方向相反時(shí)所有檢驗(yàn)方法的效能均降低，Burden檢驗(yàn)、SKAT-O和MiST降低最明顯，而SKAT、LRT和ReLRT降低較少，GenRF降低最少，LRT和ReLRT降低程度不如SKAT(表5、6)；且正負(fù)向效應(yīng)比例相差越小，效能損失越大；當(dāng)效應(yīng)方向一致時(shí)LRT和ReLRT不如其他方法，但是當(dāng)效應(yīng)方向異質(zhì)時(shí)它們優(yōu)于其他方法，效能損失相對(duì)較低。

表2 所有關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為正時(shí)的統(tǒng)計(jì)效能

表3 30%關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為負(fù)時(shí)的統(tǒng)計(jì)效能

表4 50%關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為負(fù)時(shí)的統(tǒng)計(jì)效能

表5 30%關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為負(fù)時(shí)的統(tǒng)計(jì)效能相對(duì)所有關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為正時(shí)的效能損失

表6 50%關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為負(fù)時(shí)的統(tǒng)計(jì)效能相對(duì)所有關(guān)聯(lián)位點(diǎn)效應(yīng)方向都為正時(shí)的效能損失

2.2 模擬2：可變基因片段條件下的Ⅰ型錯(cuò)誤控制和統(tǒng)計(jì)效能分析 選擇10%到90%的基因片段(對(duì)應(yīng)人類基因組不同染色體基因上可能包含不同個(gè)數(shù)的罕見變異)，樣本量固定為1 000時(shí)Ⅰ型錯(cuò)誤控制情況見表7。LRT和ReLRT在罕見變異位點(diǎn)數(shù)較少時(shí)能夠很好控制Ⅰ型錯(cuò)誤，但當(dāng)位點(diǎn)數(shù)逐漸增大時(shí)，會(huì)產(chǎn)生比名義檢驗(yàn)水準(zhǔn)更高的Ⅰ型錯(cuò)誤；而其他方法基本都能較好控制Ⅰ型錯(cuò)誤。

表7 可變基因片段和固定樣本量條件下的Ⅰ型錯(cuò)誤控制

K：包含在不同基因片段區(qū)域內(nèi)的平均罕見變異個(gè)數(shù)。

選擇10%到90%的基因片段，將其中30%的位點(diǎn)設(shè)為關(guān)聯(lián)位點(diǎn)，樣本量為1 000，假設(shè)50%的位點(diǎn)效應(yīng)b=-0.4|lgMAF|，而剩下50%的位點(diǎn)效應(yīng)b=0.4|lgMAF|。結(jié)果顯示，隨著罕見位點(diǎn)數(shù)的增加，每種檢驗(yàn)方法的統(tǒng)計(jì)效能也相應(yīng)增加，但在相同位點(diǎn)數(shù)情況下，因?yàn)榧僭O(shè)正向和負(fù)向效應(yīng)同時(shí)存在，因此LRT和ReLRT的效能最高，其次是MiST和SKAT，而Burden檢驗(yàn)的效能最低(表8)。

表8 可變基因片段條件下的統(tǒng)計(jì)效能

K：包含在該亞片段區(qū)域內(nèi)的平均罕見變異位點(diǎn)數(shù)。

2.3 運(yùn)算時(shí)間 假設(shè)檢驗(yàn)1 000個(gè)基因，罕見變異位點(diǎn)數(shù)分別為10、30、50、80和100，樣本量為2 000。計(jì)算在2.40 GHz和2.00 Gb內(nèi)存?zhèn)€人電腦上的執(zhí)行時(shí)間。表9顯示，Burden檢驗(yàn)的計(jì)算速度最快，SKAT、SKAT-O和MiST次之，LRT和ReLRT的速度較之更慢，GenRF最慢。

表9 運(yùn)算時(shí)間 s

3 討論

作者將病例-對(duì)照設(shè)計(jì)的罕見變異關(guān)聯(lián)性分析轉(zhuǎn)化為在廣義混合效應(yīng)模型框架下的似然比方差成分檢驗(yàn)。雖然這種轉(zhuǎn)化在概念上簡單，但是存在很多難題，包括參數(shù)估計(jì)、對(duì)數(shù)似然值的計(jì)算和無效分布的推導(dǎo)。作者通過PQL算法估計(jì)參數(shù)，然后在工作應(yīng)變量的基礎(chǔ)上建立偽數(shù)據(jù)；將偽數(shù)據(jù)看作原始數(shù)據(jù)并且將工作應(yīng)變量看作連續(xù)表型，這樣問題就轉(zhuǎn)化為了線性混合效應(yīng)模型的方差成分檢驗(yàn)。模擬結(jié)果表明，SKAT和MiST等方法在小樣本時(shí)相對(duì)保守，而似然比檢驗(yàn)?zāi)軌蚝芎玫乜刂脾裥湾e(cuò)誤，且具有更高的統(tǒng)計(jì)效能，但運(yùn)算速度較慢，其計(jì)算主要消耗在估計(jì)logistic混合效應(yīng)模型和模擬無效分布。

似然比檢驗(yàn)依賴于PQL算法，而PQL算法事實(shí)上是一種有偏估計(jì)，尤其是在稀疏logistic回歸中[12-15]。在稀疏情況下工作應(yīng)變量的正態(tài)近似效果差，導(dǎo)致參數(shù)估計(jì)偏倚。但從模擬結(jié)果來看這種有偏性對(duì)Ⅰ型錯(cuò)誤影響不大，尤其是在罕見變異位點(diǎn)數(shù)較少時(shí)。然而，罕見變異在位點(diǎn)數(shù)比較多時(shí)偽似然比檢驗(yàn)會(huì)變得更加自由，產(chǎn)生比名義檢驗(yàn)水準(zhǔn)更高的Ⅰ型錯(cuò)誤。其原因可能是：①點(diǎn)數(shù)增加意味著模型中隨機(jī)效應(yīng)的維數(shù)過大，此時(shí)PQL算法可能變得更加有偏。②構(gòu)建的偽數(shù)據(jù)對(duì)正態(tài)分布的近似性較差。因此，有必要進(jìn)一步研究似然比檢驗(yàn)在參數(shù)無偏估計(jì)時(shí)的表現(xiàn)，以及提高其計(jì)算速度。

[1]MANOLIO A,COLLINS S,COX J,et al.Finding the missing heritability of complex diseases[J].Nature,2009,461(7265):747

[2]EICHLER EE,FLINT J,GIBSON G,et al.Missing heritability and strategies for finding the underlying causes of complex disease[J].Nat Rev Genet,2010,11(6):446

[3]MAHER B.Personal genomes: The case of the missing heritability[J].Nature,2008,456(7218):18

[4]BODMER W,BONILLA C.Common and rare variants in multifactorial susceptibility to common diseases[J].Nat Genet,2008,40(6):695

[5]GIBSON G.Rare and common variants: twenty arguments[J].Nat Rev Genet,2012,13(2):135

[6]DICKSON P,WANG K,KRANTZ I,et al.Rare variants create synthetic genome-wide associations[J].PLoS Biol,2010,8(1):e1000294

[7]VENABLES WN,RIPLEY BD.Modern applied statistics with S[M].4th ed.Heidelberg:Springer-Verlag,2002.

[8]SCHAFFNER SF,FOO C,GABRIEL S,et al.Calibrating a coalescent simulation of human genome sequence variation[J].Genome Res,2005,15(11):1576

[9]WU MC,LEE S,CAI T,et al.Rare-variant association testing for sequencing data with the sequence kernel association test[J].Am J Hum Genet,2011,89(1):82

[10]DUCHESNE P,DE MICHEAUX PL.Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods[J].Comput Stat Data Anal,2010,54(4):858

[11]HE ZH,ZHANG M,ZHAN XW,et al.Modeling and testing for joint association using a genetic random field model[J].Biometrics,2014,70(3):471

[12]BRESLOW NE,CLAYTON DG.Approximate inference in generalized linear mixed models[J].J Am Stat Assoc,1993,88(421):9

[13]LIN X,BRESLOW NE.Bias correction in generalized linear mixed models with multiple components of dispersion[J].J Am Stat Assoc,1996,91(435):1007

[14]LIN X.Variance component testing in generalised linear models with random effects[J].Biometrika,1997,84(2):309

[15]WANG N,LIN XH,GUTIERREZ RG,et al.Bias analysis and SIMEX approach in generalized linear mixed measurement error models[J].J Am Stat Assoc,1998,93(441):249

(2016-05-17收稿 責(zé)任編輯王 曼)

Development and simulations of likelihood ratio test for rare variants association analysis in case control studies

ZENGPing1),ZHAOYang2),CHENFeng2),JINYingliang1),ZHANGYaodong1)

1)DepartmentofEpidemiologyandBiostatistics,SchoolofPublicHealth,XuzhouMedicalUniversity,Xuzhou，Jiangsu221004 2)DepartmentofBiostatistics,SchoolofPublicHealth,NanjingMedicalUniversity,Nanjing210009

rare variant;association study;logistic mixed model;likelihood ratio test

Aim: To develop a likelihood ratio test for rare variant association analysis in case control studies. Methods: The likelihood ratio test was constructed under the framework of logistic mixed models; a new pseudo-data set was obtained via the working response in the PQL algorithm, which transformed the problem of testing rare variants association in case control studies into the problem of testing rare variants association with continuous under the framework of linear mixed models. The Monte Carlo simulation was conducted to evaluate the proposed test and to compare with existing methods.Results: The simulations showed that all the methods could control the type Ⅰ error correctly. When the effects of the causal rare variants were in the same direction, the Burden test, SKAT-O and MiST were most powerful; while if both positive and negative effects were present, the likelihood ratio test outperformed the others.Conclusion: The likelihood ratio test developed using the logistic mixed models and the PQL algorithm can be applicable to rare variant association analysis in case control studies.

10.13705/j.issn.1671-6825.2017.02.004

*國家自然科學(xué)基金項(xiàng)目 81402765；國家統(tǒng)計(jì)局全國統(tǒng)計(jì)科學(xué)研究項(xiàng)目 2014LY112

R195.1