SKAT與懲罰回歸模型兩階段策略在基因組關聯研究中的應用*

廣東藥科大學公共衛生學院流行病與衛生統計學系(510310)

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

?

SKAT與懲罰回歸模型兩階段策略在基因組關聯研究中的應用*

廣東藥科大學公共衛生學院流行病與衛生統計學系(510310)

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

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

SKAT 懲罰回歸模型 基因組關聯研究

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

本研究基于同時在組水平和位點水平上推斷的思路，嘗試將SKAT與LASSO和EN聯合，應用兩階段策略進行關聯性分析，并與單水平及成組懲罰模型(cMCP和GEL)進行比較以評價各類方法的性能，為遺傳關聯研究方法學選擇提供依據和指導。

統計方法原理

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

1.SKAT

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

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

(1)

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

2.LASSO和EN

(2)

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

(3)

可看出，式(3)中當α=1時為L2懲罰，即嶺估計，當α=0時為L1懲罰，則為LASSO；實際應用中一般α取0.5。λ為調整參數，意義同式(2)[11]。研究顯示EN比LASSO可得到更加穩定、精準的預測，適用于基因微陣列等高維且存在共線性的小樣本[12]。

3.cMCP和GEL

類似SKAT集合SNPs的思想，有學者提出運用懲罰模型分析SNPs效應時也應考慮組(SNPs集)的選擇，即不僅選擇重要的成組，同時選擇組內重要的變量[13]。Breheny等在2009年提出分層懲罰的結構[6]：

(4)

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

cMCP在兩水平均使用MCP懲罰函數進行篩選，具有無偏性、稀疏性和連續性等性質，其懲罰函數和導函數分別定義為[6,14]：

(5)

式(5)中λ是決定懲罰大小的正則化參數，a是影響懲罰函數應用范圍的調節參數，如式(4)中組內和組外調節參數分別為a和b；當結局變量和協變量標準化時，推薦使用a=3[6]。式(4)中的b是組外懲罰的調節參數，為了使組水平的懲罰達到自身的最大值，可設為Qjγa/2。

在分層懲罰結構的基礎上，Breheny進一步嘗試非凸的指數懲罰函數[7]，即：

(6)

將式(6)運用于分層懲罰的結構，如式(4)中稱為group exponential lasso。可證明當τ趨于0時，式(6)退化為L1懲罰。模擬研究顯示，當組外使用指數懲罰函數(經驗值τ=1/3)，組內使用L1懲罰時(稱GEL)，該法運用變量的分組信息，同時在組水平和變量水平進行選擇，其估計準確性優于成組LASSO和cMCP[7]。

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

模擬實例分析

1.數據來源

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

2.評價方法及策略

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

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

結 果

GEL在50個數據集中模型均無法收斂，SKAT法只在基因水平上評價，因此后續在基因水平上評價8種方法，在SNPs水平上評價7種方法。

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

*50個模擬數據集分析結果；**真關聯和非關聯基因數分別為35和1106個。

各方法在基因水平的評價結果見表1。可以看到，SKAT法的平均靈敏度最高，為59.49%；SKAT+LASSO的特異度值最高，為90.60%。SKAT法的平均約登指數最高，為11.15%，其次是SKAT+EN，為8.58%。

表2 在基因水平上各方法基因選入率與基因效應的關系

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

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

*:50個模擬數據集分析結果；**:真關聯和非關聯SNPs數分別為119和531971個。

各方法在SNPs水平的評價結果見表3。可知各法靈敏度都較低，特異度均較高。相較之下，EN法與EN+SKAT的平均靈敏度最高。在特異度指標上，LASSO與其兩階段的分析策略最高。約登指數最高的為EN+SKAT與EN。

討 論

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

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

本研究還顯示除SKAT外，各方法中關聯基因的選入率均與基因效應有關。MAP4作為效應最大的基因，內含最多的真關聯SNPs；除cMCP法外，MAP4在50次試驗里幾乎均能被選出，系數值前列的SNPs亦是如此。此結果與以往運用GAW18數據評價遺傳統計方法的研究結果基本一致[19]。

限于GAW18中非獨立個體樣本量過低，本文選用了有親緣關系的研究對象，進一步研究中可在模型中納入隨機效應或采用邊際模型以解釋家庭成員表型相關。此外，從進化角度而言，稀有變異比常見變異更可能具有較強的生物學功能及遺傳效應。如定義稀有變異的閾值并施加一定權重，可能會提高統計分析的效能。

[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]張秀秀,王慧,田雙雙,等.高維數據回歸分析中基于LASSO的自變量選擇.中國衛生統計,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]曾平,趙楊,陳峰.新一代測序數據的罕見遺傳變異關聯性統計方法.中國衛生統計,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]張俊國,劉麗,李麗霞,等.懲罰廣義線性模型在遺傳關聯研究中的應用及R軟件實現.中國衛生統計,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]勾建偉.懲罰回歸方法的研究及其在后全基因關聯研究中的應用.南京醫科大學,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.

(責任編輯：郭海強)

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

國家自然科學基金(81302493)；廣東省科技廳社會發展基金(2014A020212307)；廣東省自然科學基金(2016A030313809)

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