SKAT與懲罰回歸模型兩階段策略在基因組關(guān)聯(lián)研究中的應(yīng)用*

廣東藥科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(510310)

張俊國(guó) 林志豐 劉 麗 李麗霞 楊 翌 郜艷暉△

?

SKAT與懲罰回歸模型兩階段策略在基因組關(guān)聯(lián)研究中的應(yīng)用*

廣東藥科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(510310)

張俊國(guó) 林志豐 劉 麗 李麗霞 楊 翌 郜艷暉△

目的 本研究提出兩階段分析策略，將SKAT與懲罰回歸模型聯(lián)合應(yīng)用，為遺傳關(guān)聯(lián)研究提供方法學(xué)選擇的依據(jù)和指導(dǎo)。方法 本研究使用遺傳分析工作組18的數(shù)據(jù)，分別采用SKAT，LASSO，EN，cMCP，Gel以及兩階段統(tǒng)計(jì)分析策略(SKAT+EN，SKAT+LASSO，EN+SKAT，LASSO+SKAT)進(jìn)行關(guān)聯(lián)性分析。結(jié)果 在基因水平，SKAT法的平均靈敏度與約登指數(shù)最高。除SKAT法外，其余統(tǒng)計(jì)策略的關(guān)聯(lián)基因選出率均與對(duì)結(jié)局方差解釋的比例和基因中包含SNPs的數(shù)目存在關(guān)聯(lián)。在SNPs水平，EN法與EN+SKAT的靈敏度與約登指數(shù)最高。不同的統(tǒng)計(jì)策略均能把對(duì)結(jié)局效應(yīng)貢獻(xiàn)最大的真關(guān)聯(lián)基因MAP4與SNPs篩選出來(lái)。結(jié)論 SKAT和EN聯(lián)合分析策略能夠在數(shù)百萬(wàn)SNPs中篩選主要的疾病關(guān)聯(lián)基因與SNPs，并在基因水平上統(tǒng)計(jì)推斷，有著較高靈敏度，同時(shí)還能控制嚴(yán)重的假陽(yáng)性錯(cuò)誤，為遺傳關(guān)聯(lián)研究提供了一種較為高效的統(tǒng)計(jì)分析策略。

SKAT 懲罰回歸模型 基因組關(guān)聯(lián)研究

在基因組關(guān)聯(lián)研究中，單位點(diǎn)關(guān)聯(lián)性檢驗(yàn)受限于嚴(yán)苛的多重校正顯著性水平，其效能極低[1]?；趥€(gè)體間遺傳相似性的方差分量檢驗(yàn)SKAT(sequence kernel association test)將多個(gè)SNPs(single nucleotide polymorphisms)聚集成組，不僅增加效能，還可處理位點(diǎn)間連鎖不平衡及位點(diǎn)效應(yīng)方向不同等問(wèn)題[2]。但SKAT只在組水平上進(jìn)行推斷，無(wú)法得到單個(gè)位點(diǎn)的效應(yīng)，限制了后續(xù)功能學(xué)研究的線索。由于基因組關(guān)聯(lián)數(shù)據(jù)呈現(xiàn)高維、噪音大、連鎖不平衡等特征[3]，在傳統(tǒng)最小二乘與似然估計(jì)基礎(chǔ)上引入懲罰函數(shù)是解決此類問(wèn)題的有效工具，自1996年Tibshirani提出LASSO(least absolute shrinkage and selection operator)后，懲罰回歸模型廣受關(guān)注，基于懲罰思想發(fā)展了許多統(tǒng)計(jì)學(xué)方法[4]，如LASSO和嶺回歸結(jié)合的EN(elastic net)[5]，與LASSO都可在大幅降維的同時(shí)估計(jì)單個(gè)變異的效應(yīng)?？紤]到位點(diǎn)間的相互作用和增加效能，也可同時(shí)在組水平和位點(diǎn)水平上懲罰，如cMCP(composite minimax concave penalty)[6]和GEL(group exponential lasso)[7]等，但各類方法應(yīng)用效果仍有待于進(jìn)一步研究。

本研究基于同時(shí)在組水平和位點(diǎn)水平上推斷的思路，嘗試將SKAT與LASSO和EN聯(lián)合，應(yīng)用兩階段策略進(jìn)行關(guān)聯(lián)性分析，并與單水平及成組懲罰模型(cMCP和GEL)進(jìn)行比較以評(píng)價(jià)各類方法的性能，為遺傳關(guān)聯(lián)研究方法學(xué)選擇提供依據(jù)和指導(dǎo)。

統(tǒng)計(jì)方法原理

假設(shè)有n個(gè)觀測(cè)，第i個(gè)研究對(duì)象的P個(gè)SNPs基因型數(shù)據(jù)表示為xi=(xi1,xi2,…,xiP)T,i=1,2,…,n；其中xip=0,1,2(p=1,2,…,P)分別對(duì)應(yīng)主要等位基因的純合子、雜合子以及最小等位基因的純合子，根據(jù)生物學(xué)先驗(yàn)，將待分析的P個(gè)SNPs分成J組(如以基因?yàn)榻M單位)，Qj為第j組中的SNPs個(gè)數(shù)(q=1,2,…,Qj)；進(jìn)一步假設(shè)xip已被中心化。每個(gè)研究對(duì)象有K個(gè)人口學(xué)、環(huán)境或其他混雜因素，用Zi=(zi1,zi2,…ziK)T來(lái)表示。第i個(gè)研究對(duì)象的表型狀態(tài)為yi∈R。

1.SKAT

SKAT以SNPs集(如基因、通路或ROI)為分析單位，在線性混合效應(yīng)模型的框架下，通過(guò)核函數(shù)(kernel function)量化個(gè)體間的遺傳相似性，并基于得分函數(shù)(score function)進(jìn)行方差成份檢驗(yàn)[8]，當(dāng)表型性狀為連續(xù)型變量時(shí)，模型為

yi=β0+β1zi1+…+βKziK+h(xi1,xi2,…xiQj)+εi

(1)

式(1)中β0是截距項(xiàng)，β1,…,βK表示環(huán)境或人口學(xué)特征等協(xié)變量的回歸系數(shù)，εi為隨機(jī)誤差，服從N(0，σ2)。核函數(shù)h(·)綜合了集合中所有SNPs的遺傳信息，選擇不同形式的核函數(shù)可擬合集合內(nèi)SNPs與表型的線性或非線性關(guān)聯(lián)[9-10]。

2.LASSO和EN

(2)

式(2)中λ≥0，為懲罰參數(shù)，用于控制壓縮程度。

(3)

可看出，式(3)中當(dāng)α=1時(shí)為L(zhǎng)2懲罰，即嶺估計(jì)，當(dāng)α=0時(shí)為L(zhǎng)1懲罰，則為L(zhǎng)ASSO；實(shí)際應(yīng)用中一般α取0.5。λ為調(diào)整參數(shù)，意義同式(2)[11]。研究顯示EN比LASSO可得到更加穩(wěn)定、精準(zhǔn)的預(yù)測(cè)，適用于基因微陣列等高維且存在共線性的小樣本[12]。

3.cMCP和GEL

類似SKAT集合SNPs的思想，有學(xué)者提出運(yùn)用懲罰模型分析SNPs效應(yīng)時(shí)也應(yīng)考慮組(SNPs集)的選擇，即不僅選擇重要的成組，同時(shí)選擇組內(nèi)重要的變量[13]。Breheny等在2009年提出分層懲罰的結(jié)構(gòu)[6]：

(4)

式(4)中fλ,b與fλ,a分別代表組水平與SNPs水平的懲罰函數(shù)。

cMCP在兩水平均使用MCP懲罰函數(shù)進(jìn)行篩選，具有無(wú)偏性、稀疏性和連續(xù)性等性質(zhì)，其懲罰函數(shù)和導(dǎo)函數(shù)分別定義為[6,14]：

(5)

式(5)中λ是決定懲罰大小的正則化參數(shù)，a是影響懲罰函數(shù)應(yīng)用范圍的調(diào)節(jié)參數(shù)，如式(4)中組內(nèi)和組外調(diào)節(jié)參數(shù)分別為a和b；當(dāng)結(jié)局變量和協(xié)變量標(biāo)準(zhǔn)化時(shí)，推薦使用a=3[6]。式(4)中的b是組外懲罰的調(diào)節(jié)參數(shù)，為了使組水平的懲罰達(dá)到自身的最大值，可設(shè)為Qjγa/2。

在分層懲罰結(jié)構(gòu)的基礎(chǔ)上，Breheny進(jìn)一步嘗試非凸的指數(shù)懲罰函數(shù)[7]，即：

(6)

將式(6)運(yùn)用于分層懲罰的結(jié)構(gòu)，如式(4)中稱為group exponential lasso。可證明當(dāng)τ趨于0時(shí)，式(6)退化為L(zhǎng)1懲罰。模擬研究顯示，當(dāng)組外使用指數(shù)懲罰函數(shù)(經(jīng)驗(yàn)值τ=1/3)，組內(nèi)使用L1懲罰時(shí)(稱GEL)，該法運(yùn)用變量的分組信息，同時(shí)在組水平和變量水平進(jìn)行選擇，其估計(jì)準(zhǔn)確性優(yōu)于成組LASSO和cMCP[7]。

上述懲罰回歸模型的擬合均使用坐標(biāo)下降法(coordinate descent)[15]或由其改進(jìn)的局部近似坐標(biāo)下降法(locally approximated coordinate descent)[6]。調(diào)整參數(shù)的確定則通過(guò)K折交叉驗(yàn)證、廣義交叉驗(yàn)證、無(wú)偏估計(jì)的風(fēng)險(xiǎn)分析以及BIC準(zhǔn)則等[16]。所有方法均可在R3.2.2軟件實(shí)現(xiàn)，分別調(diào)用軟件包SKAT(SKAT)glmnet(LASSO和EN)，grpreg(cMCP和GEL)。

模擬實(shí)例分析

1.數(shù)據(jù)來(lái)源

本研究數(shù)據(jù)源于遺傳分析工作組18(genetic analysis workshop 18，GAW18)[17]，是一個(gè)國(guó)際上公開(kāi)的用于研究稀有變異關(guān)聯(lián)方法的模擬數(shù)據(jù)平臺(tái)。本文選擇性別、年齡和血壓均無(wú)缺失的849個(gè)存在親緣關(guān)系的個(gè)體作為研究對(duì)象，將3個(gè)時(shí)間點(diǎn)的舒張壓(DBP)均值作為結(jié)局變量，選取對(duì)DBP方差解釋比例最大(7.79%)的3號(hào)染色體中SNPs作為自變量。先通過(guò)UCSC基因?yàn)g覽器(http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz/)對(duì)3號(hào)染色體的全部SNPs(共1215399個(gè))進(jìn)行基因標(biāo)記，基因的范圍為最小的轉(zhuǎn)錄起始位點(diǎn)到最大的轉(zhuǎn)錄結(jié)束位點(diǎn)之間；再刪除基因間與存在缺失的SNPs后，最終共標(biāo)記1141個(gè)基因的532092個(gè)SNPs，其中MAF低于1%的稀有變異占51%。納入分析的532,092個(gè)SNPs中與DBP存在真關(guān)聯(lián)的基因(SNPs)有35(119)個(gè)，真關(guān)聯(lián)位點(diǎn)解釋DBP的方差比例為7.27%。其中，對(duì)DBP效應(yīng)最大的基因?yàn)镸AP4，解釋的方差比例為6.48%。

2.評(píng)價(jià)方法及策略

GAW18提供了200個(gè)和遺傳變異有關(guān)聯(lián)的模擬表型數(shù)據(jù)集，考慮到計(jì)算負(fù)擔(dān)，本研究分析其中50個(gè)。對(duì)每個(gè)數(shù)據(jù)集，分別采用SKAT，LASSO，EN，cMCP，GEL以及兩階段聯(lián)合分析策略(SKAT+EN，SKAT+LASSO，EN+SKAT，LASSO+SKAT)進(jìn)行關(guān)聯(lián)性分析，其中聯(lián)合策略如SKAT+EN表示先用SKAT篩選，將結(jié)果陽(yáng)性的基因(SNPs)再納入EN進(jìn)行分析。由于SKAT只在基因水平上評(píng)價(jià)，故本研究中SKAT分析時(shí)選入的基因及其SNPs均視為陽(yáng)性。LASSO和EN只在SNPs水平上評(píng)價(jià)，則模型中系數(shù)不為0的SNPs及所屬基因均視為陽(yáng)性。

為評(píng)價(jià)各統(tǒng)計(jì)方法識(shí)別真關(guān)聯(lián)基因(SNPs)以及剔除非關(guān)聯(lián)基因的能力；本研究在基因水平和SNPs水平上分別計(jì)算50個(gè)數(shù)據(jù)集中各方法的平均靈敏度(%)、特異度(%)和約登指數(shù)(%)。其中靈敏度定義為：各方法篩選出的真關(guān)聯(lián)基因(SNP)的數(shù)目除以基因組中實(shí)際關(guān)聯(lián)基因(SNP)的總數(shù)；特異度定義為：各方法剔除的非關(guān)聯(lián)基因(SNP)的數(shù)目除以基因組中實(shí)際非關(guān)聯(lián)基因(SNP)的總數(shù)；約登指數(shù)定義為靈敏度+特異度-1。進(jìn)一步采用Spearman相關(guān)系數(shù)(rs)評(píng)價(jià)各方法選入的基因中基因解釋DBP方差的比例、基因內(nèi)包含的真關(guān)聯(lián)SNPs數(shù)目與基因選入率(50次運(yùn)算中某基因選入的比例)的關(guān)系。其中SKAT法、cMCP和GEL均以基因作為組水平。因?yàn)锽onferroni校正過(guò)于苛刻，故SKAT法采用FDR法進(jìn)行多重校正。懲罰回歸中調(diào)整參數(shù)確定均采用10折交叉驗(yàn)證。所有統(tǒng)計(jì)分析中均忽略有親緣關(guān)系個(gè)體間的家庭相關(guān)。

結(jié) 果

GEL在50個(gè)數(shù)據(jù)集中模型均無(wú)法收斂，SKAT法只在基因水平上評(píng)價(jià)，因此后續(xù)在基因水平上評(píng)價(jià)8種方法，在SNPs水平上評(píng)價(jià)7種方法。

表1 在基因水平上各類方法的平均靈敏度(%)、特異度(%)和約登指數(shù)(%)*

*50個(gè)模擬數(shù)據(jù)集分析結(jié)果；**真關(guān)聯(lián)和非關(guān)聯(lián)基因數(shù)分別為35和1106個(gè)。

各方法在基因水平的評(píng)價(jià)結(jié)果見(jiàn)表1?？梢钥吹剑琒KAT法的平均靈敏度最高，為59.49%；SKAT+LASSO的特異度值最高，為90.60%。SKAT法的平均約登指數(shù)最高，為11.15%，其次是SKAT+EN，為8.58%。

表2 在基因水平上各方法基因選入率與基因效應(yīng)的關(guān)系

如表2所示，除SKAT法外，其余方法基因選入率均與基因效應(yīng)有關(guān)，基因解釋DBP方差比例越大、基因內(nèi)真關(guān)聯(lián)SNPs數(shù)越多，則越容易被選入。對(duì)DBP效應(yīng)最大的基因MAP4，除cMCP外，其他方法在50次模擬中選入率均較高(98%～100%)。

表3 在SNPs水平上各類方法的平均靈敏度(%)、特異度(%)和約登指數(shù)(%)*

*:50個(gè)模擬數(shù)據(jù)集分析結(jié)果；**:真關(guān)聯(lián)和非關(guān)聯(lián)SNPs數(shù)分別為119和531971個(gè)。

各方法在SNPs水平的評(píng)價(jià)結(jié)果見(jiàn)表3。可知各法靈敏度都較低，特異度均較高。相較之下，EN法與EN+SKAT的平均靈敏度最高。在特異度指標(biāo)上，LASSO與其兩階段的分析策略最高。約登指數(shù)最高的為EN+SKAT與EN。

討 論

在全基因組關(guān)聯(lián)研究中，從浩瀚如煙的遺傳變異中篩選與疾病存在關(guān)聯(lián)的少數(shù)病因變異對(duì)統(tǒng)計(jì)方法提出巨大的挑戰(zhàn)。SKAT與懲罰回歸模型是近年來(lái)熱門的遺傳統(tǒng)計(jì)方法；SKAT理論上有吸引力且計(jì)算快捷，在基因水平上靈敏度高，可篩選更多的真關(guān)聯(lián)基因。懲罰回歸模型具有良好的預(yù)測(cè)精度與穩(wěn)定性，能夠在數(shù)十萬(wàn)SNPs中大幅度壓縮，挑選出與疾病關(guān)聯(lián)性最強(qiáng)的基因與SNPs；但與以往研究[18]一致，單獨(dú)使用懲罰回歸模型仍會(huì)產(chǎn)生大量的假陽(yáng)性，需要在獨(dú)立人群中進(jìn)一步驗(yàn)證。本研究將兩法結(jié)合，結(jié)果表明，可將更多真關(guān)聯(lián)基因(SNPs)納入后續(xù)分析(SKAT+LASSO、SKAT+EN)或?qū)Y(jié)果做進(jìn)一步篩選(LASSO+SKAT、EN+SKAT)。在基因水平上，SKAT的性能最高，其余方法中兩階段策略的性能略優(yōu)于單一策略。在SNPs水平上，EN法與EN+SKAT的性能略高于其它策略。

盡管兩水平懲罰模型理論合理，但本研究顯示在基因組關(guān)聯(lián)研究中，cMCP的性能指標(biāo)均低于其余統(tǒng)計(jì)分析策略，并未表現(xiàn)出優(yōu)勢(shì)。該法在SNPs水平上的性能與LASSO相差無(wú)幾，對(duì)變量系數(shù)的兩次壓縮并不能很大程度上減少方差，反而引入不必要的偏差，損失了更多的真關(guān)聯(lián)SNPs[7]。此外，GEL法在本研究所有數(shù)據(jù)中模型均無(wú)法收斂，當(dāng)變量數(shù)遠(yuǎn)遠(yuǎn)大于觀測(cè)數(shù)的時(shí)候，該法可能無(wú)法在較低的β值上實(shí)現(xiàn)模型擬合，此時(shí)模型無(wú)法識(shí)別或接近于奇異，系數(shù)的路徑也不存在。

本研究還顯示除SKAT外，各方法中關(guān)聯(lián)基因的選入率均與基因效應(yīng)有關(guān)。MAP4作為效應(yīng)最大的基因，內(nèi)含最多的真關(guān)聯(lián)SNPs；除cMCP法外，MAP4在50次試驗(yàn)里幾乎均能被選出，系數(shù)值前列的SNPs亦是如此。此結(jié)果與以往運(yùn)用GAW18數(shù)據(jù)評(píng)價(jià)遺傳統(tǒng)計(jì)方法的研究結(jié)果基本一致[19]。

限于GAW18中非獨(dú)立個(gè)體樣本量過(guò)低，本文選用了有親緣關(guān)系的研究對(duì)象，進(jìn)一步研究中可在模型中納入隨機(jī)效應(yīng)或采用邊際模型以解釋家庭成員表型相關(guān)。此外，從進(jìn)化角度而言，稀有變異比常見(jiàn)變異更可能具有較強(qiáng)的生物學(xué)功能及遺傳效應(yīng)。如定義稀有變異的閾值并施加一定權(quán)重，可能會(huì)提高統(tǒng)計(jì)分析的效能。

[1]Gang P,Li L,Hoicheong S,et al.Gene and pathway-based second-wave analysis of genome-wide association studies.European Journal of Human Genetics Ejhg,2010,18(1):111-117.

[2]Wu M,Lee S,Cai T,et al.Rare-variant association testing for sequencing data with the sequence kernel association test.American Journal of Human Genetics,2011,89(1):82-93.

[3]張秀秀,王慧,田雙雙,等.高維數(shù)據(jù)回歸分析中基于LASSO的自變量選擇.中國(guó)衛(wèi)生統(tǒng)計(jì),2013,30(6):922-926.

[4]Tibshirani R.Regression Shrinkage and Selection via the Lasso.Journal of the Royal Statistical Society,1996,58(1):267-288.

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

[6]Breheny P,Huang J.Penalized methods for bi-level variable selection.Statistics & Its Interface,2009,2(3):369-380.

[7]Breheny P.The group exponential lasso for bi-level variable selection.Biometrics,2015,71(3):731-740.

[8]曾平,趙楊,陳峰.新一代測(cè)序數(shù)據(jù)的罕見(jiàn)遺傳變異關(guān)聯(lián)性統(tǒng)計(jì)方法.中國(guó)衛(wèi)生統(tǒng)計(jì),2015,32(6):1091-1096.

[9]Brown MP,Grundy WN,Lin D,et al.Knowledge-based analysis of microarray gene expression data by using support vector machines.Proceedings of the National Academy of Sciences.2000,97(1):262-267.

[10]Liu D,Ghosh D,Lin X.Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernel machine regression via logistic mixed models.Bmc Bioinformatics.2008,9(14):292.

[11]張俊國(guó),劉麗,李麗霞,等.懲罰廣義線性模型在遺傳關(guān)聯(lián)研究中的應(yīng)用及R軟件實(shí)現(xiàn).中國(guó)衛(wèi)生統(tǒng)計(jì),2016，33(4):582-586.

[12]Hesterberg T,Choi NH,Meier L,et al.Least angle and1 penalized regression:A review.Statistics Surveys.2008,2.

[13]Huang J,Ma S,Xie H,et al.A group bridge approach for variable selection.Biometrika,2009,96(2):339-355.

[14]Zhang CH.Nearly unbiased variable selection under minimax concave penalty.Annals of Statistics,2010,38(2):894-942.

[15]Friedman J,Hastie T,Tibshirani R.Regularization paths for generalized linear models via coordinate descent.Journal of Statistical Software,2009,33(1):1-22.

[16]Fu WJ.Nonlinear GCV and quasi-GCV for shrinkage models.Journal of Statistical Planning & Inference,2005,131(2):333-347.

[17]Laura A,Dyer TD,Peralta JM,et al.Data for Genetic Analysis Workshop 18:human whole genome sequence,blood pressure,and simulated phenotypes in extended pedigrees.Bmc Proceedings,2014,8(1):1-9.

[18]勾建偉.懲罰回歸方法的研究及其在后全基因關(guān)聯(lián)研究中的應(yīng)用.南京醫(yī)科大學(xué),2014.

[19]Cordell HJ.Summary of Results and Discussions From the Gene-Based Tests Group at Genetic Analysis Workshop 18.Genetic Epidemiology,2014,38 Suppl 1(S1):S44-S48.

(責(zé)任編輯：郭海強(qiáng))

Two-steps Strategies Jointing SKAT with Penalized Regression and their Application in Genome-wide Association Study

Zhang Junguo,Lin Zhifeng,Liu Li,et al

(DepartmentofEpidemiologyandBiostatistics,SchoolofPublicHealth,GuangdongPharmaceuticalUniversity(510310),Guangzhou)

Objective This study proposes two-stage analysis strategy to combine the advantages of two types of methods in order to provide a method guidance for the genetic association study.Methods SKAT,LASSO,EN and two-stage strategies(SKAT+EN，SKAT+LASSO，EN+SKAT，LASSO+SKAT)as well as bi-level variable selection models(cMCP,Gel)are used in the data of the genetic analysis workshop 18 to compare their application performance.Results At the gene level show that the method of SKAT has the highest average sensitivity and average Youden index.The rate of gene of these statistical methods except the method of SKAT are associated with the number of SNPs within the gene and the proportion of explained variance of DBP.The result at the SNP level indicate that the method of EN has highest sensitivity.The highest Youden index is counted by EN+SKAT method and the second is EN method.The gene of MAP4 and SNPs that is the largest contribution to DBP all selected by the various statistical analysis.Conclusion The combination of the methods of EN and SKAT could screen few number variants that associate with phenotypes in big data.This methods not only has high sensitivity but also has restraint false positives,it could provide some clues for the future studies of genetic mechanisms.

SKAT;Penalized regression;Genome-wide association study

國(guó)家自然科學(xué)基金(81302493)；廣東省科技廳社會(huì)發(fā)展基金(2014A020212307)；廣東省自然科學(xué)基金(2016A030313809)

△通信作者：郜艷暉，E-mail:gao_yanhui@163.com。