鄧舒方，王昱泉，胡躍清

(復旦大學 生命科學學院 生物統(tǒng)計學與計算生物學系，上海 200433)

得益于高通量測序技術的快速發(fā)展，我們獲得了大量全基因組范圍內(nèi)的遺傳變異數(shù)據(jù)，而如何有效地利用它們進行基因與疾病的關聯(lián)分析是生物統(tǒng)計學的重要任務.在早期的全基因組關聯(lián)分析(Genome-Wide Association Study, GWAS)研究中，人們一般以單個單核苷酸多態(tài)性(Single Nucleotide Polymorphism, SNP)為單位進行疾病的關聯(lián)分析[1-2]，而隨著研究的深入，人們發(fā)現(xiàn)復雜疾病的發(fā)生往往涉及多個位點，其中一些位點的效應可能較弱，此時，單SNP分析可能會漏篩一些致病位點[3-7].為了提高檢驗功效，需要對多個SNP進行整合分析： 如果存在一個致病位點，由于相鄰位點間存在連鎖不平衡，故相鄰位點也會與疾病關聯(lián)，將這一段區(qū)域內(nèi)的位點進行聯(lián)合分析能最大程度地提高檢驗功效.基于這一思想，Wang等[8]提出了一種針對病例對照數(shù)據(jù)的多位點聯(lián)合關聯(lián)分析方法SLIDE(the test Statistic incorporating Linkage Disequilibrium)，該方法能夠有效利用位點間的連鎖不平衡(Linkage Disequilibrium, LD)信息對統(tǒng)計量的方差協(xié)方差矩陣進行計算，從而提高檢驗功效.

然而，遺傳學研究數(shù)據(jù)中常常包含一些協(xié)變量信息，如年齡[9]、種族[10-11]等，這些協(xié)變量信息在進行遺傳位點與性狀的關聯(lián)性檢測中有著重要作用，如果不對這些協(xié)變量信息進行適當處理，可能會帶來分析結果的偏差或是檢驗功效的降低[12].當協(xié)變量既與預測變量相關也與響應變量相關時，若在探究預測變量和響應變量的關聯(lián)性過程中不考慮因此帶來的混淆影響，會增加分析結果的假陽性率[13].此外，調(diào)整協(xié)變量在某種程度上能夠通過降低統(tǒng)計量方差來增加預測模型的準確度[14-15].因此，在遺傳關聯(lián)分析中，對協(xié)變量進行調(diào)整是有意義的.

在多位點關聯(lián)分析中，傳統(tǒng)的調(diào)整協(xié)變量的方法主要有兩類： 第一類是將響應變量對預測變量及協(xié)變量同時回歸[16]，如SKAT[17]及SKAT-O[18]對應的調(diào)整協(xié)變量方法;第二類是將響應變量對協(xié)變量進行回歸得到殘差，將殘差再與預測變量進行關聯(lián)分析[19]，如傳統(tǒng)的殘差回歸[20-21]以及SSU/aSPU[22]對應的調(diào)整協(xié)變量方法.這兩類方法都依賴于假設的疾病產(chǎn)生模型，普適性不高.故我們基于SLIDE這一不依賴于模型的多位點非參關聯(lián)檢測方法，提出一種對基因型用傾向得分進行逆概率加權[23-24]的調(diào)整協(xié)變量關聯(lián)分析方法SLIDEa.

為了檢驗SLIDEa方法的表現(xiàn)，我們計算了不同隨機模擬場景下該方法與文獻中現(xiàn)有方法的第一類錯誤率以及檢驗功效.模擬結果表明，它在基因與協(xié)變量交互作用存在或基因間效應方向相反時均具有較好表現(xiàn).隨后我們將SLIDEa應用至酗酒相關遺傳學合作研究數(shù)據(jù)集COGA中，找到了線粒體外膜脂代謝酶基因OPA3、催產(chǎn)素受體基因OXTR等數(shù)個與酒精成癮性相關的基因[25-27]，這表明了SLIDEa方法的有效性并能為后續(xù)的生物學研究提供理論指引.

1 符號和方法

1.1 現(xiàn)有方法介紹

SKAT方法是Wu等[17]于2011年提出的一種關聯(lián)分析方法，它基于病例對照數(shù)據(jù)的如下模型

其中：G是預測變量；D是二分類響應變量；Z是協(xié)變量.原假設下所有預測變量與響應變量無關聯(lián)，其統(tǒng)計量為

SKAT-O方法是Lee等[18]于2012年提出的，它其實是SKAT的統(tǒng)計量與Burden檢驗統(tǒng)計量的加權和.

Qρ=(1-ρ)QSKAT+ρQBurden，

其中：ρ為合并方法Burden所占權重，0≤ρ≤1;Pρ為統(tǒng)計量Qρ的P值.在原假設下，Qoptimal服從混合卡方分布.考慮到協(xié)變量Z容易納入到回歸模型中，我們用adj_SKAT/adj_SKAT-O分別表示調(diào)整了協(xié)變量Z的方法.

aSPU是Pan等[22]于2009年提出的一種關聯(lián)分析方法.對于病例對照研究，aSPU從模型

得到得分統(tǒng)計量向量

其中：ξj為每個位點得分統(tǒng)計量的權重，1≤j≤m.而aSPU統(tǒng)計量為

其中：γ為權重；PSPU(γ)為檢驗統(tǒng)計量TSPU(γ)的P值；Γ是γ的一個取值范圍，如在本文的隨機模擬中，Γ={1，2，4，∞}.特別地，當Γ={2}時，aSPU就是SSU.

在協(xié)變量存在的情況下，adj_aSPU與adj_SSU首先將響應變量對協(xié)變量做回歸，得到殘差項，隨后將殘差項當成新的響應變量得到調(diào)整協(xié)變量后的得分統(tǒng)計量向量，后續(xù)與aSPU/SSU計算方式相同.adj_resi與adj_aSPU類似，先將響應變量對協(xié)變量做回歸，得到殘差項，隨后將殘差項當成新的響應變量對預測變量做回歸.

1.2 多位點關聯(lián)分析方法SLIDE介紹

對病例對照組數(shù)據(jù)，Wang等[8]在回溯性框架下提出了一個能同時對m個遺傳變異位點進行關聯(lián)分析的統(tǒng)計量TSLIDE，其U向量及其方差協(xié)方差矩陣表達式如下：

其中：nA表示患病個體數(shù)；nU為健康個體數(shù)；第i個個體的患病狀態(tài)記為Di，其中Di=1表示患病，Di=0表示健康；Gi=(Gi1，Gi2，…，Gim)T為第i個個體的m個位點基因型，其中Gij表示第i個個體在第j個位點的基因型，Gij={0，1，2}，它表示較小等位基因在該位點上的個數(shù).U統(tǒng)計量反映了多位點基因型在病例組和對照組之間的差異.在原假設H0： 這m個位點均與疾病無關聯(lián)成立的條件下，患病人群與健康人群的多位點基因型均值無差異，從而E(U)=0，對應的方差協(xié)方差矩陣為

Cov(G)=(σjk)m×m，σjj=2pj(1-pj)，σjk=2δjk， 1≤j≠k≤m，

其中：G表示研究所在的總體中某個體在m個位點的基因型；pj為第j個位點的較小等位基因頻率，1≤j≤m;δjk為第j、k個位點間的LD值，δjk=pjk-pjpk，pjk為在一條染色單體上第j、k位點基因均為較小頻率等位基因的概率，檢驗統(tǒng)計量為

在原假設下，TSLIDE漸近服從自由度為m的卡方分布，從而方便得到它的P值.

1.3 基于SLIDE的逆概率加權調(diào)整協(xié)變量方法

非對角元(j，k)(j≠k)為

2 隨機模擬研究

2.1 數(shù)據(jù)產(chǎn)生模型

我們先通過模型產(chǎn)生10萬個個體的協(xié)變量、基因型信息和性狀值，再依據(jù)0/1型性狀從相應人群中各隨機抽500個個體作為病例和對照用于下一步的多位點關聯(lián)分析.數(shù)據(jù)產(chǎn)生模型具體分為3步：

(i) 產(chǎn)生p維協(xié)變量信息Z.對于第i個個體，1≤i≤n=100 000，Zi=(Zi1，Zi2，…，Zip)T，其中Zi1，Zi2，…，Zip相互獨立，均服從Bernoulli(0.5)分布.

(ii) 由Zi通過有序logistic回歸產(chǎn)生第i個個體基因型Gi=(Gi1，Gi2，…，Gim)T，其中位點j的基因型Gij生成模型是

其中： 1≤j≤m，m取值為5或50；(ε1，ε2，…，εm)T來自于均值為0的m維正態(tài)分布；εj、εk間的協(xié)方差Cov(εj，εk)=ρ|j-k|，1≤j≠k≤m；ρ的大小反映了位點間的相關程度；模型中的μgj根據(jù)對應位點的最小等位基因頻率(Minor Allele Frequency, MAF)及樣本協(xié)變量信息計算得到，即μgj需要滿足等式E(Gij)=2×MAFj，MAFj表示第j個位點較小等位基因在人群中的頻率；αj=(αj1，αj2，…，αjp)T表示p維協(xié)變量對第j個基因位點的影響.

(iii) 由Zi、Gi通過logistic回歸產(chǎn)生性狀Di：

其中：β0反映人群發(fā)病率；βg為位點對性狀的m×1效應向量；βz為協(xié)變量對性狀的p×1效應向量；βgz為位點與協(xié)變量的交互項對性狀的m×p效應矩陣.

對于第一類錯誤率，我們共考慮SNP個數(shù)(5或50)、協(xié)變量維數(shù)p(1或2)和位點類型(全為常見變異(Common Variant, CV)或CV與罕見變異(Rare Variant, RV)交替出現(xiàn))的6種不同組合情形.CV的MAF從0.1至0.4中隨機產(chǎn)生;RV的MAF從0.005至0.01中隨機產(chǎn)生;每種情形中又根據(jù)協(xié)變量對性狀是否有效應、ρ=0或ρ=0.5分為4種子情形;每個情形下重復模擬次數(shù)為1 000，顯著性水平為0.05;為了計算第一類錯誤率，設β0=ln0.1、位點對性狀無效應、位點與協(xié)變量交互項對性狀也無效性，詳細參數(shù)設置見表1，其中：α=(α1，α2，…，αm)為p×m參數(shù)矩陣.

表1 原假設下的參數(shù)設置Tab.1 Parameter setting under H0

隨后，我們根據(jù)SNP個數(shù)、協(xié)變量維數(shù)、位點MAF、位點與協(xié)變量交互作用方向的不同設置了8種情形來比較各種方法的功效： 當SNP個數(shù)為5時，假設前兩個位點為致病位點；當SNP個數(shù)為50時，假設第1、2、6、7個位點為致病位點，其中常見致病位點對性狀效應值的取值范圍為0.05～0.3、罕見變異位點效應值范圍為0.09～1.09;情形1～情形4中設置了5個SNP及一維協(xié)變量，考慮了是否僅含常見變異、位點與協(xié)變量交互作用方向相同/相反下不同組合的情形，用以探究不同位點類型、不同基因與協(xié)變量相互作用效應方向下各種方法的表現(xiàn);情形5～情形6將協(xié)變量擴展為二維，用以探究在二維協(xié)變量情況下各方法功效差異;情形7～情形8將SNP個數(shù)擴大為50，用以探究各方法在多SNP個數(shù)下的表現(xiàn);在情形1～情形8中，又考慮了致病位點對性狀效應方向相同/相反兩種子情形，如在5個SNP中，βg=(ln1.05，ln1.05，0，0，0)T是當位點效應方向相同的情況，那么與之對應βg=(ln1.05，ln(1/1.05)，0，0，0)T則為位點效應方向相反情況;具體參數(shù)設置見表2，模擬研究重復次數(shù)為1 000，顯著性水平為0.05.

表2 計算功效的參數(shù)設置Tab.2 Parameter setting for power evaluation

2.2 隨機模擬結果

表3展示了各方法及其對應的調(diào)整協(xié)變量方法在原假設情形1～情形6下的第一類錯誤率.從中不難發(fā)現(xiàn)當βz不為零時，未進行協(xié)變量調(diào)整的方法的第一類錯誤率明顯大于0.05，這說明當協(xié)變量既與基因型又與性狀相關時，不考慮協(xié)變量直接進行基因與性狀的關聯(lián)性分析會控制不住第一類錯誤率，故在這種情形下，對協(xié)變量進行調(diào)整是減小假陽性率的有效手段.且當SNP個數(shù)增加至50個時，SLIDEa方法仍能控制住第一類錯誤率.

表3 原假設中各情形下各方法的第一類錯誤率Tab.3 Type I error rate for each scenario under H0

接著在α、βz均不為零，即協(xié)變量既與基因型有關也與性狀有關的情形下對各調(diào)整協(xié)變量方法進行功效比較.圖1(見第774頁)和圖2(見第775頁)為對立假設下8種情形各方法的功效比較.當SNP個數(shù)為5、協(xié)變量為一維時，如情形1～情形4中結果所示，無論不同位點對性狀效應方向一致還是相反，在位點與協(xié)變量存在交互作用的情況下，本文提出的SLIDEa方法均相較于其它多位點調(diào)整協(xié)變量關聯(lián)分析方法具有優(yōu)越性，并在不同位點與協(xié)變量交互項效應方向相反時它仍能保持明顯優(yōu)勢;當關聯(lián)分析位點包含罕見變異位點時，SLIDEa方法雖然仍能保持優(yōu)勢，但這種優(yōu)勢在關聯(lián)分析位點全為常見變異位點時表現(xiàn)得更為明顯;當增加協(xié)變量維數(shù)至二維、保持SNP個數(shù)不變時，如情形5～情形6中結果所示，當位點與協(xié)變量存在交互作用的情況下，SLIDEa的表現(xiàn)均優(yōu)于其他方法;當SNP個數(shù)增加至50個、位點全為常見變異時，無論協(xié)變量維數(shù)是一維還是二維，SLIDEa同樣具有較好表現(xiàn).

圖1 不同致病位點對性狀效應方向相同時對立假設情形1～情形8下各方法功效比較Fig.1 Power comparison for scenarios 1 to 8 under Ha while the effect directions of these loci are identical注： 該表展示了對立假設情形下當不同致病位點對性狀效應方向相同時，情形1～情形8下各個方法的功效;其中βg為第一個位點對性狀的回歸系數(shù).

圖2 不同致病位點對性狀效應方向相反時對立假設情形1～情形8下各方法功效比較Fig.2 Power comparison for scenarios 1 to 8 under Ha while the effect directions of these loci are opposite注： 該表展示了對立假設情形下當不同致病位點對性狀效應方向相反時，情形1～情形8下各個方法的功效;其中βg為第一個位點對性狀的回歸系數(shù).

綜上可得，SLIDEa適用于基因和環(huán)境存在交互作用時的基因與性狀間關聯(lián)分析，無論位點對性狀的效應方向相同或相反、不同位點與協(xié)變量交互項效應方向相同或相反，它均能保持其優(yōu)勢.TSLIDEa統(tǒng)計量在構建過程中未依賴任何疾病發(fā)生模型，而其他方法，如adj_resi、adj_ aSPU等雖然也是基于無模型假設，但在調(diào)整協(xié)變量時，需要將性狀對協(xié)變量進行回歸，從而得到性狀的殘差，再用殘差與基因進行關聯(lián)分析.因此，在用其它方法調(diào)整協(xié)變量時，若是模型假定與真實模型不一致，會導致統(tǒng)計量功效降低.SLIDEa在調(diào)整協(xié)變量和關聯(lián)分析中都不依賴于實際疾病發(fā)生模型，因此該方法是較為穩(wěn)健的，受疾病真實發(fā)病模型的影響較小.

3 實際數(shù)據(jù)分析

酗酒相關遺傳學合作研究數(shù)據(jù)(The Collaborative Study on the Genetics of Alcoholism, COGA)的參與者來自于美國7個不同地區(qū)，它收集了詳細的酒精依賴診斷所需信息.我們分析的是COGA中一部分依據(jù)《精神障礙診斷與統(tǒng)計手冊第4版》(DSM-Ⅳ)標準對進行酒精成癮性評估劃分為病例組及對照組后篩選出來的病例對照研究數(shù)據(jù)，包含共1 945人的全基因組測序數(shù)據(jù)和年齡、性別、每日吸煙量等30個潛在酗酒相關協(xié)變量因子.除去測序缺失的數(shù)據(jù)，我們最終得到1 913個個體共包含107萬個遺傳位點、18 946個基因的測序數(shù)據(jù).接下來，我們在NCBI(National Center for Biotechnology Information)上篩選出304個已報導的與酒精成癮性相關的基因，用SLIDEa、aSPU、SKAT和SKAT-O對應的協(xié)變量調(diào)整方法對其中在COGA測序數(shù)據(jù)中出現(xiàn)的267個基因逐一進行關聯(lián)分析.我們將年齡和家庭酗酒史作為協(xié)變量進行調(diào)整，來探究在給定年齡及家族史的情況下與酗酒相關的關聯(lián)性位點.

表4展示了用SLIDEa得到的前10位顯著基因和對應P值及其他方法得到的P值.從該表中可以看出，各個方法均能發(fā)現(xiàn)這些基因，其P值也在相近水平，但是SLIDEa的P值會稍低一些.通過文獻檢索，我們可以得出這些基因與酗酒成癮性具有一定關聯(lián).例如，其中的CHRNB2為神經(jīng)元煙堿類受體，能夠與乙酰膽堿或尼古丁結合控制鈉、鉀離子通過細胞膜，這一基因與尼古丁依賴高度相關[28]，后續(xù)全基因組關聯(lián)分析中也發(fā)現(xiàn)它還與酒精成癮性相關聯(lián)[29-30];OPA3基因編碼蛋白參與線粒體內(nèi)膜凈化[31]，而線粒體內(nèi)膜參與酒精代謝[32]，已發(fā)表的GWAS研究也有發(fā)現(xiàn)OPA3與酒精成癮性相關[25];而OXTR作為G蛋白偶聯(lián)受體，也在多個基因與酒精成癮性GWAS研究中被檢測到[26-27].這些基因與酒精成癮性的關系也進一步說明了SLIDEa的可靠性.

表4 COGA數(shù)據(jù)利用SLIDEa得到的前10位顯著基因及所有方法的P值Tab.4 P-values for top 10 significant genes obtained from SLIDEa on COGA data

這10個基因由Metascape工具(網(wǎng)址： http：∥metascape.org/)得到的注釋結果詳見表5.將這10個基因輸入PANTHER網(wǎng)站(網(wǎng)址： http：∥www.pantherdb.org/)做基因富集分析，可以得到4個GO通路，分別是GO： 0042220(P=1.24×10-2)，GO： 0051952(P=4.54×10-4)，GO： 0043269(P=4.61×10-2)，GO： 0003008(P=2.61×10-2).這4個基因通路分別與可卡因、胺轉(zhuǎn)運、離子轉(zhuǎn)運和系統(tǒng)過程相關，和酒精成癮性的關系尚不明確，有待進一步研究.

表5 COGA數(shù)據(jù)利用SLIDEa得到的前10位顯著基因注釋結果Tab.5 Annotation result of top 10 significant genes obtained from SLIDEa on COGA data

4 總結與討論

已有的研究結果表明了在協(xié)變量既與基因有關也與性狀有關時調(diào)整協(xié)變量的必要性.為了消除關聯(lián)分析時協(xié)變量帶來的混淆影響、降低關聯(lián)分析方法的假陽性率，本文在多位點聯(lián)合關聯(lián)分析方法SLIDE的基礎上，提出了一種對基因型用傾向得分進行逆概率加權的調(diào)整協(xié)變量關聯(lián)分析方法SLIDEa，其中的協(xié)變量可以是多維的.模擬結果顯示，當致病位點為常見變異時，與aSPU、SKAT等調(diào)整協(xié)變量方法相比，在基因?qū)π誀钚喾椿虼嬖诨蚺c協(xié)變量的交互作用時，SLIDEa具有一定優(yōu)勢;而當致病位點有較多罕見變異時，它的效果不是很理想，這是因為在等位基因頻率幾乎為0時，用樣本數(shù)據(jù)估計U的方差會出現(xiàn)問題，從而使SLIDEa的功效降低.SLIDEa的另一優(yōu)勢在于它的構建過程不依賴于疾病產(chǎn)生模型，這使得它的適用范圍較廣.隨后，我們將提出的方法應用于酗酒遺傳學實際數(shù)據(jù)COGA的分析中，找到的顯著性基因可以用來做進一步的生物學基礎研究.

本方法仍有幾個方面值得改進： 1) 需要提升它在存在較多罕見變異位點時的效果，與SLIDE相比，它無法利用外部大樣本數(shù)據(jù)來估算統(tǒng)計量的方差，所以在致病位點全為罕見變異的情況下，它的表現(xiàn)不是很理想.2) 當關聯(lián)分析SNP位點數(shù)目m大、樣本量不大時，存在數(shù)據(jù)稀疏性，此時統(tǒng)計量可能會不服從卡方分布.此外，統(tǒng)計量方差中涉及到與m2同階多的條件期望的計算;故當關聯(lián)分析位點數(shù)目過多時，可以根據(jù)位點間的位置將多位點劃分成數(shù)目較小的一些片段，得到這些片段的P值后整合得到整個基因的P值.