曹紅艷 曾 平 李 治 崔躍華，4張巖波△

懲罰廣義估計(jì)方程在縱向數(shù)據(jù)基因關(guān)聯(lián)分析中的應(yīng)用*

曹紅艷1曾 平2李 治3崔躍華1，4張巖波1△

目的 探討懲罰廣義估計(jì)方程（pGEE）在縱向數(shù)據(jù)基因關(guān)聯(lián)分析的應(yīng)用，為縱向數(shù)據(jù)基因關(guān)聯(lián)分析提供方法學(xué)參考。方法 以小鼠糖尿病發(fā)病相關(guān)的數(shù)量性狀位點(diǎn)識別為例，分別采用廣義估計(jì)方程（GEE）和pGEE進(jìn)行分析。結(jié)果 pGEE篩選出糖尿病發(fā)病關(guān)聯(lián)位點(diǎn)，為分子生物學(xué)研究提供了重要的候選位點(diǎn)。結(jié)論 pGEE能有效的實(shí)現(xiàn)高維縱向數(shù)據(jù)的變量選擇，識別出有意義的關(guān)聯(lián)位點(diǎn)。

懲罰廣義估計(jì)方程 縱向數(shù)據(jù) SCAD 基因關(guān)聯(lián)分析 數(shù)量性狀位點(diǎn)

隨著分子生物學(xué)測序技術(shù)的不斷發(fā)展，基因關(guān)聯(lián)分析已經(jīng)成為了復(fù)雜疾病研究的最重要手段之一，成功的識別了多種人類復(fù)雜疾?。ㄈ纰蛐吞悄虿?，高血壓等）的遺傳變異位點(diǎn)［1-2］。為了更好地了解復(fù)雜疾病和遺傳變異之間的關(guān)系，越來越多的基因關(guān)聯(lián)分析通過縱向的方式進(jìn)行，對同一觀察對象的某觀察指標(biāo)在不同時間上進(jìn)行重復(fù)測量，數(shù)據(jù)形式表現(xiàn)為非獨(dú)立結(jié)構(gòu)?？v向數(shù)據(jù)基因關(guān)聯(lián)分析相比于橫斷面研究，可以研究復(fù)雜性狀隨時間變化的關(guān)系，從而增強(qiáng)統(tǒng)計(jì)效能，提高遺傳變異對復(fù)雜疾病的解釋程度［3-5］。

縱向基因數(shù)據(jù)和一般的縱向數(shù)據(jù)相比，更為復(fù)雜，高維度，不同SNP之間存在復(fù)雜的LD結(jié)構(gòu)，表現(xiàn)為強(qiáng)相關(guān)性和多維共線性；同時，大部分的SNP為冗余信息，只有少量的SNP為關(guān)聯(lián)位點(diǎn)且信號強(qiáng)度弱，將導(dǎo)致參數(shù)估計(jì)和統(tǒng)計(jì)推斷的準(zhǔn)確性和有效性大大降低，高維縱向基因關(guān)聯(lián)分析面臨著巨大的挑戰(zhàn)。因此，發(fā)展新的變量選擇方法尤為重要。近年來，基于懲罰的變量選擇方法備受關(guān)注，如Lasso［6］，SCAD［7］，自適應(yīng)Lasso［8］等。懲罰方法能有效地應(yīng)用于高維數(shù)據(jù)分析，通過收縮將弱效應(yīng)估計(jì)為0，參數(shù)估計(jì)和變量選擇同時進(jìn)行，廣泛的應(yīng)用于不同模型的變量選擇。經(jīng)典的縱向數(shù)據(jù)分析中，廣義估計(jì)方程（generalized estimating equations，GEE）［9-11］只需定義偽得分方程（quasi -score equations），當(dāng)作業(yè)相關(guān)矩陣設(shè)置不正確時，仍然能得到一致性的參數(shù)估計(jì)，在應(yīng)用中獨(dú)占優(yōu)勢。因此，Wang，Zhou和Qu發(fā)展了基于SCAD的懲罰廣義估計(jì)方程（penalized generalized estimating equations，pGEE）［12］，不僅保持了GEE的重要特性，同時將GEE推廣到高維數(shù)據(jù)分析，適用于協(xié)變量個數(shù)p隨樣本例數(shù)n同階變化的情況，即P＝O（n）。本文以小鼠的糖尿病發(fā)病相關(guān)的數(shù)量性狀位點(diǎn)（quantitative trait locus，QTL）識別為例，進(jìn)行pGEE分析，識別出小鼠的糖尿病發(fā)病相關(guān)的QTL，為QTL分析提供方法支持。

資料與方法

1．資料來源

數(shù)據(jù)來源于2型糖尿病發(fā)病相關(guān)的QTL遺傳研究［13］，由于肥胖是2型糖尿病主要危險因素，Reifsnyder等將肥胖且有糖尿病傾向的NZO（new zealand obese）／HILt小鼠與瘦且無糖尿病傾向的NON（nonobese nondiabetic）／Lt小鼠進(jìn)行遠(yuǎn)交（outcross），產(chǎn)生的F1代再與親本NON／Lt小鼠回交（backcross）生成多基因?qū)嶒?yàn)小鼠模型。對203例雄性實(shí)驗(yàn)小鼠，測量其在4、8、12、16、20、24周的體重，同時，在已知QTL位點(diǎn)周圍20 cM范圍內(nèi)的高頻標(biāo)記中，挑選了83個微衛(wèi)星位點(diǎn)，研究微衛(wèi)星位點(diǎn)對體重的影響，數(shù)據(jù)形式如表1所示。

表1 體重和微衛(wèi)星標(biāo)記數(shù)據(jù)集形式

2．方法

（1）廣義估計(jì)方程

在縱向數(shù)據(jù)基因關(guān)聯(lián)分析中，設(shè)Yit為個體i在時間t的觀察值，Xit為p維協(xié)變量向量，令Yi＝（Yi1，…Yit）T為個體i的觀察向量，Xi＝（Xi1，…，Xit）T為個體i的協(xié)變量矩陣，則縱向數(shù)據(jù)分析的模型框架可表達(dá)如下：

E（Yit）＝μi，g（μi）＝X′iβ

基于全似然函數(shù)的參數(shù)估計(jì)方法，只有在一些特殊情況下，如當(dāng)Yit服從多元正態(tài)分布時，才能直接算出，對于其他的分布，全似然函數(shù)的計(jì)算非常復(fù)雜。GEE采用類似于廣義線性模型的得分方程的方法，定義偽得分方程為：

（2）懲罰廣義估計(jì)方程

基于懲罰的變量選擇方法有很多，F(xiàn)an和Li指出一個好的懲罰函數(shù)估計(jì)值應(yīng)具備無偏性、稀疏性和連續(xù)性，即Oracle性質(zhì)［7］。SCAD懲罰能保留較大的系數(shù)，同時將較小的系數(shù)收縮為0，具有Oracle性質(zhì)，其懲罰函數(shù)pλ（θ）的導(dǎo)數(shù)如下：

其中，I為指示函數(shù)，a為預(yù)先選擇的常數(shù)，F(xiàn)an和Li推薦a＝3．7。

對GEE的得分函數(shù)進(jìn)行SCAD懲罰，得到pGEE的懲罰表達(dá)式為：

采用minorization-maximization算法得到（2）式的漸進(jìn)表達(dá)式，

其中Sj（β＾）為S（β＾）的第j個得分函數(shù)，ε的取值較小，一般取ε＝10-6。進(jìn)一步采用New ton-Raphson算法對（3）進(jìn)行迭代，計(jì)算得到pGEE的參數(shù)估計(jì)值。pGEE估計(jì)值依賴于懲罰參數(shù)λ，在此，采用交叉驗(yàn)證（cross validation，CV）選擇最優(yōu)λ。

（3）統(tǒng)計(jì)分析采用R軟件中的“PGEE 1．4”軟件包，可實(shí)現(xiàn)GEE和pGEE分析，GEE分析也可采用“geepack”包?！癙GEE 1．4”包含3個函數(shù)：MGEE，CVfit和PGEE?！癕GEE”函數(shù)用于GEE分析，“CVfit”用于計(jì)算交叉驗(yàn)證的最優(yōu)懲罰參數(shù)λ，得到最優(yōu)λ后，采用“PGEE”函數(shù)計(jì)算pGEE的參數(shù)估計(jì)值。適用于測量結(jié)果呈正態(tài)分布、二項(xiàng)分布、Poisson分布以及Gamma分布的情況，為縱向數(shù)據(jù)基因關(guān)聯(lián)分析提供了強(qiáng)有力的分析工具。

本例中，以體重為應(yīng)變量，以83個微衛(wèi)星位點(diǎn)、測量時間及時間的平方為協(xié)變量，共85個協(xié)變量，X＝［1，time，time2，x1，x2，…，x83］，其中xj編碼為0和1，這是典型的縱向數(shù)據(jù)基因關(guān)聯(lián)分析。對所有變量均進(jìn)行標(biāo)準(zhǔn)化，采用正態(tài)鏈接函數(shù)的GEE和pGEE進(jìn)行分析。由于每只實(shí)驗(yàn)小鼠在等距的時間間隔內(nèi)重復(fù)測量了6次，兩次測量相隔時間越長，相應(yīng)的相關(guān)關(guān)系越小，因此，假定作業(yè)相關(guān)矩陣為1階自回歸結(jié)構(gòu)。截距項(xiàng)不進(jìn)行懲罰。通過“CVfit”函數(shù)選擇出最優(yōu)的懲罰參數(shù)為λ＝0．0212，進(jìn)一步利用“PGEE”函數(shù)進(jìn)行擬合，得到pGEE參數(shù)估計(jì)結(jié)果。GEE分析在“geepack”軟件包中實(shí)現(xiàn)，GEE和pGEE參數(shù)設(shè)置均采用軟件中的默認(rèn)設(shè)置。

結(jié) 果

1．實(shí)驗(yàn)小鼠體重一般情況

將每只小鼠的體重測量值按時間進(jìn)行連接，得體重隨時間變化趨勢圖，如圖1?？梢娦∈笤?4周內(nèi)的個體生長趨勢總體上一致，呈先快后慢的趨勢，但也存在個體間的差異，個體增長幅度可能不一致，其差別隨時間的變化而增加。進(jìn)一步分析每個微衛(wèi)星位點(diǎn)的單獨(dú)效應(yīng)，對每個位點(diǎn)分別擬合GEE，得到單個位點(diǎn)的效應(yīng)系數(shù)，如圖2，單個位點(diǎn)系數(shù)分布在［-0．176，0．099］之間，大部分位點(diǎn)的效應(yīng)非常弱。

圖1 小鼠不同周次體重變化趨勢圖

圖2 單個QTL系數(shù)圖

2．GEE和pGEE分析

在考慮了時間效應(yīng)的情況下，觀察GEE和pGEE對83個微衛(wèi)星位點(diǎn)的參數(shù)估計(jì)散點(diǎn)圖（圖3），pGEE將大部分具有微弱效應(yīng)的位點(diǎn)收縮為0，有效的進(jìn)行了降維，而未懲罰的GEE卻不滿足稀疏性。GEE采用Wald卡方檢驗(yàn)進(jìn)行統(tǒng)計(jì)推斷，選出了14個有意義的QTL（檢驗(yàn)水準(zhǔn)），pGEE篩選出11個QTL，其參數(shù)估計(jì)結(jié)果見表2。比較兩方法共同篩選出的QTL參數(shù)估計(jì)值，可看出pGEE估計(jì)值明顯小于GEE，從QTL效應(yīng)的微弱性而言，pGEE估計(jì)出的QTL效應(yīng)更符合分子生物學(xué)特點(diǎn)。

圖3 83個微衛(wèi)星位點(diǎn)的GEE和pGEE參數(shù)估計(jì)散點(diǎn)圖

在pGEE篩選出的11個QTL中，根據(jù)每個QTL系數(shù)的正負(fù)，認(rèn)為D1M it211，D5M it158，D6M it275及D15M it29的突變?yōu)槲ｋU因素，將導(dǎo)致小鼠肥胖，其余7個QTL的變異對小鼠肥胖具有保護(hù)效應(yīng)。

肥胖與2型糖尿病的發(fā)病顯著相關(guān)，將誘發(fā)糖尿病的發(fā)生。肥胖的發(fā)生、發(fā)展受多個微效基因及其復(fù)雜的交互作用控制，同時也受環(huán)境因子的調(diào)控，是一種復(fù)雜的多因素、多基因疾?。?4］。本例中通過pGEE分析篩選出的QTL，和Reifsnyder采用單因素方差分析和置換檢驗(yàn)識別的體重相關(guān)的7個第一染色體位點(diǎn)相比［13］，共同識別的位點(diǎn)為D1M it211，D1M it411，D1M it76。陳峰采用線性混合效應(yīng)模型（linear mixed model，LMM）的Lasso和SCAD懲罰對小鼠體重QTL進(jìn)行了識別［15］，分別識別出14和3個QTL，和本文結(jié)果相比較，基于Lasso的LMM、基于SCAD的LMM以及本文的pGEE，三種不同懲罰模型共同識別出了D1M it411和D5M it158位點(diǎn)，提示這兩個位點(diǎn)的重要性，同時pGEE識別的其他9個QTL，也將為糖尿病發(fā)病相關(guān)QTL分子生物學(xué)研究提供重要的候選位點(diǎn)。

表2 GEE和pGEE參數(shù)估計(jì)值

總之，pGEE成功的對高維縱向基因數(shù)據(jù)進(jìn)行了降維，其參數(shù)估計(jì)具有oracle性質(zhì)，有效的實(shí)現(xiàn)了高維縱向基因數(shù)據(jù)的變量選擇，篩選出了有意義的遺傳變異。

討 論

pGEE分析借助作業(yè)相關(guān)矩陣考慮了不同時間點(diǎn)測量值之間的內(nèi)部相關(guān)關(guān)系，通過基于SCAD的懲罰方法，實(shí)現(xiàn)了高維縱向數(shù)據(jù)的變量選擇。本文以糖尿病縱向數(shù)據(jù)的QTL分析為例，識別出糖尿病發(fā)病相關(guān)的QTL，說明了pGEE在縱向數(shù)據(jù)基因關(guān)聯(lián)性分析中應(yīng)用的科學(xué)性和可行性。

基于Wald卡方檢驗(yàn)的GEE統(tǒng)計(jì)推斷，其前提是在保證其他變量不變的情況下，估計(jì)變量的偏效應(yīng)，當(dāng)維度較高且樣本量較小時，基于偏效應(yīng)的統(tǒng)計(jì)推斷無法揭示變量的真實(shí)效應(yīng)?；趹土P的方法參數(shù)估計(jì)和變量選擇同時進(jìn)行，pGEE中采用的SCAD懲罰具有一致性和oracle性質(zhì)，保證了非零系數(shù)變量的正確推斷和選擇。在高維縱向基因關(guān)聯(lián)分析中，pGEE的變量選擇明顯優(yōu)于GEE的參數(shù)估計(jì)和統(tǒng)計(jì)推斷。

pGEE變量選擇的一致性和稀疏性依賴于懲罰參數(shù)的選擇，基于交叉驗(yàn)證的懲罰參數(shù)選擇方法容易出現(xiàn)過擬合，從而將無效變量選為有意義的變量，出現(xiàn)錯誤選擇［16］。Hyunkeun等通過數(shù)據(jù)模擬和真實(shí)數(shù)據(jù)分析，指出基于貝葉斯信息準(zhǔn)則（bayesian information criterion，BIC）的懲罰參數(shù)選擇方法得到的參數(shù)估計(jì)值具有一致性［17］。因此，今后可進(jìn)一步研究pGEE框架下的BIC懲罰參數(shù)選擇方法，選擇出最優(yōu)的懲罰參數(shù)，以得到更為確切的變量選擇結(jié)果。需要注意的是，經(jīng)pGEE篩選出的有意義的位點(diǎn)，還需要結(jié)合分子生物學(xué)實(shí)驗(yàn)進(jìn)一步驗(yàn)證其生物學(xué)意義。

總之，pGEE作為縱向數(shù)據(jù)經(jīng)典方法GEE的懲罰模型，一方面，延續(xù)了GEE在縱向數(shù)據(jù)分析中的重要優(yōu)勢：只需設(shè)定一階和二階矩以及作業(yè)相關(guān)矩陣，避免了高維縱向數(shù)據(jù)分析中更為復(fù)雜的全似然函數(shù)計(jì)算；當(dāng)作業(yè)相關(guān)矩陣指定不當(dāng)時，仍能保持參數(shù)估計(jì)的一致性。另一方面，pGEE將GEE推廣到高維數(shù)據(jù)分析，既能考慮縱向數(shù)據(jù)之間的相關(guān)關(guān)系，又能有效的進(jìn)行降維，識別出高維縱向基因數(shù)據(jù)的關(guān)聯(lián)位點(diǎn)，為高維縱向數(shù)據(jù)分析提供了強(qiáng)有力的分析方法。隨著分子生物學(xué)技術(shù)成本的降低，縱向基因數(shù)據(jù)必將日益增多，pGEE將在縱向數(shù)據(jù)基因關(guān)聯(lián)性分析中發(fā)揮及其重要的作用。

［1］M ccarthy M I，Zeggini E．Genome-w ide association studies in type 2 diabetes．Curr Diab Rep，2009，9（2）：164-171．

［2］Ehret GB．Genome-w ide association studies：contribution of genom ics to understanding blood pressure and essential hypertension．Curr Hypertens Rep，2010，12（1）：17-25．

［3］Sitlani CM，Rice K M，Lum ley T，et al．Generalized estimating equations for genome-w ide association studies using longitudinal phenotype data．Stat Med，2015，34（1）：118-130．

［4］Furlotte NA，Eskin E，Eyheramendy S．Genome-w ide association mapping w ith longitudinal data．Genet Epidem iol，2012，36（5）：463-471．

［5］ShiG，Rao D．Ignoring temporal trends in genetic effects substantially reduces power of quantitative trait linkage analysis．Genetic epidem iology，2008，32（1）：61-72．

［6］Tibshirani R．Regression shrinkage and selection via the lasso．Journal of the Royal Statistical Society．Series B（Methodological），1996：267-288．

［7］Fan J，Li R．Variable selection via nonconcave penalized likelihood and its oracle properties．Journal of the American statistical Association，2001，96（456）：1348-1360．

［8］Zou H．The adaptive lasso and its oracle properties．Journal of the A-merican statistical association，2006，101（476）：1418-1429．

［9］Liang KY，Zeger SL．Longitudinal data analysis using generalized linearmodels．Biometrika，1986，73（1）：13-22．

［10］陳啟光．縱向研究中重復(fù)測量資料的廣義估計(jì)方程分析．中國衛(wèi)生統(tǒng)計(jì)，1995，12（1）：22-25．

［11］李洪艷，譚珊，高曉，等．基于廣義估計(jì)方程的嬰兒超重的影響因素分析．中國衛(wèi)生統(tǒng)計(jì)，2016，33（2）：222-225．

［12］Wang L，Zhou J，Qu A．Penalized generalized estimating equations for high-dimensional longitudinal data analysis．Biometrics，2012，68（2）：353-360．

［13］Reifsnyder PC，Churchill G，Leiter EH．Maternal environment and genotype interact to establish diabesity in m ice．Genome Res，2000，10（10）：1568-1578．

［14］Reifsnyder PC，Leiter EH．Deconstructing and reconstructing obesityinduced diabetes（diabesity）in m ice．Diabetes，2002，51（3）：825-832．

［15］陳峰．線性混合效應(yīng)模型的懲罰變量選擇．中國衛(wèi)生信息管理雜志，2014，11（3）：278-284．

［16］Wang H，Li R，Tsai CL．Tuning parameter selectors for the smoothly clipped absolute deviation method．Biometrika，2007，94（3）：553-568．

［17］Cho H，Qu A．Model selection for correlated data with diverging number of parameters．Statistica Sinica，2013，23：901-927．

（責(zé)任編輯：張 悅）

The Application of Penalized Generalized Estimating Equations in Genetic Association w ith Longitudinal Data

Cao Hongyan，Zeng Ping，Li Zhi，et al（Department of Health Statistics，ShanxiMedical University（030001），Taiyuan）

Objective To explore the application of penalized generalized estimating equations（pGEE）in genetic association w ith longitudinal data，and provide new statistical solutions for longitudinal genetic data analysis．M ethods We applied the generalized estimating equations（GEE）and pGEEmethods to a type II diabetes dataset for quantitative trait locus identification．Results Several loci associated w ith the development of type IIdiabeteswere identified using the pGEEmethod，providing important candidatemakers for further biological validation．Conclusion The pGEEmethod provides a powerful tool for high dimensional longitudinal genetic association studies．

Penalized generalized estimating equations；Longitudinal data；SCAD；Genetic association studies；Quantitative trait locus

*：國家自然科學(xué)基金資助項(xiàng)目（30972553，31371336）

1．山西醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)教研室（030001）

2．徐州醫(yī)科大學(xué)流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室

3．中北大學(xué)體育學(xué)院

4．美國密西根州立大學(xué)統(tǒng)計(jì)與概率系

△通信作者：張巖波，E-mail：sxmuzyb＠126．com