基于分而治之的Lasso方法對海量基因數(shù)據(jù)的研究與驗證

蘭曉然，張 灝,2

（1.太原理工大學 數(shù)學學院，太原 030024；2.亞利桑那大學 數(shù)學系，美國 亞利桑那州 圖森 85721）

0 引言

近來醫(yī)學發(fā)展尤其迅猛，但仍有許多遺傳疾病的病因和發(fā)病機制尚未明確，多數(shù)學者認為與遺傳因素及環(huán)境因素等多種因素共同作用有關(guān)，遺傳與環(huán)境交互作用起到重要作用，屬于多個易感基因參與的復雜疾病[1]。因此，定位與性狀或疾病相關(guān)聯(lián)的位點，對醫(yī)學研究人員了解性狀和一些疾病的遺傳機理至關(guān)重要。與此同時，在基因?qū)W研究中會出現(xiàn)高維海量數(shù)據(jù)集，如何在大量基因中選擇出強相關(guān)的位點變量，是變量選擇研究要迫切解決的問題。

基因位點的選擇研究，要求變量選擇方法應(yīng)保持的特點有：（1）可解釋性，能很好地解釋基因位點對疾病的作用。（2）強相關(guān)位點選擇的穩(wěn)定性。（3）盡量控制計算機的消耗時間。這些問題也是變量選擇至關(guān)重要的問題，但原來已有的方法如：最優(yōu)子集選擇、逐步回歸、嶺回歸等只具有其中部分性質(zhì)。因此，如何解決這些問題，成為變量選擇界一大熱點。在Frank提出ridge Regression和Bireman提出Nonnegative Garrote算法的基礎(chǔ)上，Tibshirani于1996年提出了一種新的算法即Lasso，此方法成功的應(yīng)用于COX模型。雖然Lasso方法很好地彌補了傳統(tǒng)方法上的不足[2]，但是針對海量基因數(shù)據(jù)，仍會出現(xiàn)計算機開銷大，運行時間過長的問題。

針對上述問題，本文提出一種基于分而治之的Lasso算法，再通過采用Q-Q圖對篩選結(jié)果進行驗證，并選用卡方和Fisher檢驗對篩選出的位點進行顯著性檢驗，從結(jié)果表現(xiàn)改進的Lasso方法可以有效地對海量基因數(shù)據(jù)進行變量選擇，挑選出與疾病強相關(guān)的致病位點，并極大地縮小了計算機耗費時間。

1 Lasso方法

Lasso方法是1996年Tibshirani提出的既能對變量進行選擇[3],又能得出參數(shù)估計值的一種變量選擇方法。

考慮如下普通線性模型：

其中Y=(y1，y2，…，，yn)T為響應(yīng)變量，n為樣本容量，X=(X1，X2，…，Xn)為p維預測變量，Xi=(Xi1，Xi2，…，Xin)T，i=1，2，…，n，β=(β1，β2，…，βp)T稀疏，即β1，β2，…，βp有很多系數(shù)為零，ε=(ε1，ε2，…，εn)T是正態(tài)分布隨機誤差，即ε～N(0，σ2In) 。 假 設(shè) 觀 測 數(shù) 據(jù) (yi，xij)，i=1，2，…，n，j=1，2，…，p已經(jīng)過中心標準化處理，即：

在下文中，除特別說明外，數(shù)據(jù)（X,Y）均為已經(jīng)過中心標準化處理的數(shù)據(jù)。

對固定非負數(shù)義，Lasso方法定義如下:

Lasso方法不僅繼承了傳統(tǒng)方法的優(yōu)點，而且是十分有效的變量選擇算法，使得其在統(tǒng)計學中受到更為廣泛的關(guān)注與研究。Lars算法是由Efron等提出的，目前在R語言中有Lars程序包，直接調(diào)入就能進行變量選擇的有關(guān)計算，非常方便實用。

2 分而治之方法

針對高維海量的基因數(shù)據(jù)，變量具有稀疏性[4]。用Lasso進行變量選擇，會使計算機消耗過大，計算時間過長。所以，本文提出了分而治之方法，它不僅能夠更好的排除錯誤模型選擇帶來的偽相關(guān)，而且可以極大地降低計算時間。變量選擇的時間復雜度一致于O(napb)，a＞1，b≥0[5]。

假設(shè)參數(shù)數(shù)量p相當大，定義，是稀疏的。假設(shè)總數(shù)據(jù)集的大小為n，將其劃分為K份，則第k個子集有nk個觀測值：( )xk，i，yk，i，i=1，…，nk。 記，則對第k個子集（k=1，…，K）的對數(shù)似然函數(shù)：

相應(yīng)的懲罰估計為:

其中ρ(β;λk)訓練參數(shù)λk的懲罰函數(shù)，可參見Fan和Lv（2011）[6]。

3 實驗分析

3.1 實驗數(shù)據(jù)

本文數(shù)據(jù)來自16年研究生數(shù)學建模，數(shù)據(jù)詳細描述如表1所示。

表1 屬性變換表

研究關(guān)聯(lián)基因的選擇問題已成為一個熱點話題，很多關(guān)于MS的GWAS和Meta[7]的文章陸續(xù)發(fā)表。

3.2 實驗過程

試驗中，對基因數(shù)據(jù)運用Lasso與分而治之方法進行分析比較。首先觀察兩種方法在計算時間消耗的大小，說明分而治之算法能很大的節(jié)省時間，降低消耗。然后對篩選出的變量進行驗證。具體步驟如下：

（1）時間消耗對比

表2 時間消耗對比表

由表2可以看出，經(jīng)過split-and-conquer操作后，預測精度影響不大，但明顯縮短了運行時間。以下將N/5為分塊標準在進一步研究。

（2）調(diào)整參數(shù)lambda的確定

對lambda的格點值，進行5折交叉驗證，選取交叉驗證均方誤差誤差最小的lambda值。然后，按照得到的lambda值，用全部數(shù)據(jù)重新擬合模型（見下頁圖1、圖2）。

可以看到最佳的lambda取值就是在紅色曲線的最低點處，Lasso與分而治之方法選擇出來的參數(shù)都是在誤差接近0.24的值約為-3。

（3）變量篩選

Lasso和分而治之算法采用的算法是循環(huán)坐標下降法，因為遺傳病為0～1分類問題，所以采用Logistic回歸模型。

由圖3可看出，隨著橫軸L1范數(shù)的增加來調(diào)整系數(shù)，可以看出在橫軸為15的范圍內(nèi)只有7個變量的系數(shù)不為0，位點rs2273298第一個被選擇出來。由此可以得出分而治之選擇出來7種致病位點分別為rs12036216、rs7368252、rs2273298、rs932372、rs7522344、rs12133956、rs2143810。

圖1 Lasso方法lambda選擇圖

圖2 分而治之方法lambda選擇圖

圖3 分而治之方法位點選擇圖

3.3 實驗結(jié)果分析與檢驗

3.3.1 變量選擇結(jié)果分析

為了減小隨機性帶來的誤差，本文對變量選擇進行了100次循環(huán)，則各變量在100次中被選次數(shù)如表3所示。

使用Lasso方法對數(shù)據(jù)進行篩選，在100次的篩選中，位點rs12036216、rs7368252、rs2273298、rs932372每次都會選中，而位點rs3013045、rs11573253、rs6683624選中次數(shù)較少，位點rs6683624僅被選擇了1次。相對于Lasso方法，分而治之方法選擇的位點更集中，位點rs12036216、rs7368252、rs2273298、rs932372在100次中都被選擇，其余的位點rs7522344、rs12133956、rs2143810也多次被選擇。說明分而治之方法在保持穩(wěn)定性的基礎(chǔ)上，選擇的區(qū)間更緊湊。從理論上來說，分而治之方法在對每份子數(shù)據(jù)集進行變量選擇時，去除與類屬性不相關(guān)的變量，其中包含部分冗余的變量，但是再對所有變量進行系統(tǒng)整合再進行選擇，可以刪除冗余變量。因而，可以保證變量選擇的有效性。

表3 致病位點被選次數(shù)統(tǒng)計表

3.3.2 Q-Q Plot驗證結(jié)果

Q-Q Plot主要是用來估計數(shù)量性狀觀測值與預測值之間的差異。在GWAS研究中，Q-Q Plot的X和Y軸主要是代表各個SNP的-lg(P.values)[8]。預測的線是一條從原點發(fā)出的45°角的虛線，實際觀測值則是標的實心點。在GWAS研究中，如果哪個SNP點出現(xiàn)了較大的偏離，則認為這個SNP位點的觀測值的偏離是由這個SNP突變所產(chǎn)生的遺傳作用造成的。Q-Q Plot的具體算法為：將P值由大到小進行排序，進而運用公式計算得出橫縱坐標，計算分別為公式為：驗證結(jié)果如圖4所示。

圖4 Q-Q圖

圖4中，細線表示零假設(shè)下的期望值，黑點表示對不同遺傳位點和性狀關(guān)聯(lián)性分析的結(jié)果。通過對位點觀測值的統(tǒng)計量與期望統(tǒng)計量進行比較，可以看出遺傳位點間是否存在群體分層現(xiàn)象，并驗證是否存在一些位點對遺傳性疾病的具有顯著的影響。針對Q-Q Plot通常呈現(xiàn)的結(jié)果，圖4顯示，藍線和紅線下端重合，表明遺傳位點不存在群體分層，但上端存在一些偏差，表明存在對遺傳疾病具有顯著影響的遺傳位點。并且由圖可看出，有7個遺傳性位點的觀測值統(tǒng)計量偏離期望統(tǒng)計量，表明有7個遺傳性位點對遺傳病具有顯著影響，這一結(jié)論恰好與運用分而治之分析出的7個關(guān)聯(lián)性大的位點相一致。

3.3.3 卡方檢驗和Fisher檢驗

卡方檢驗是根據(jù)卡方分布的原理[9，10]，計算檢驗統(tǒng)計量卡方的值，再與卡方分布的臨界值比較，確定P值并做出推斷，卡方值計算的基本公式為：

其中A為實際頻數(shù)，T為理論頻數(shù)。

在群體遺傳學中卡方檢驗是應(yīng)用最多的一種假設(shè)檢驗方法，常用來判斷某事物的頻數(shù)分布是否符合某一理論分布，若符合即可按該理論分布來處理。就本實驗而言，用Lasso方法和分而治之方法分析出的數(shù)據(jù)關(guān)聯(lián)性最高的遺傳位點，可以用卡方檢驗來驗證其實驗得出的數(shù)據(jù)是否具有顯著性。

在二分類問題中，由于某些類別的例數(shù)較少，在本文采用Fisher檢驗，假設(shè)有2×2列聯(lián)表，如表4所示。

表4 2×2列聯(lián)表

如果固定行和列，那么在零假設(shè)條件下出現(xiàn)在四格表中的各數(shù)值分別為n11，n12，n21，n22，假設(shè)邊緣頻數(shù)n1·，n2·，n·1，n·2和n··都是固定的，在A和B獨立或沒有齊性的零假設(shè)下，對任意的i，j，nij服從超幾何分布為：

分別用卡方和Fisher方法來檢驗位點與某遺傳病的關(guān)聯(lián)性，檢驗結(jié)果如表5所示。

表5 顯著性檢驗

通過表5分析可得，卡方檢驗和Fisher檢驗在5%的顯著性水平下，均拒絕原假設(shè)，通過5%顯著性檢驗，表明運用Lasso、分而治之方法篩選出的7個遺傳位點對遺傳位點具有顯著的影響作用。

綜合分析，本文在篩選變量時，運用Lasso、分而治之方法篩選出的位點，再運用Q-Q圖、卡方以及Fisher方法進行檢驗。結(jié)果表明，分而治之方法能夠有效地應(yīng)用于海量的基因數(shù)據(jù)中，對其進行變量選擇，不僅節(jié)省了計算時間，較小計算機消耗，而且保證了選擇變量與類標簽的強相關(guān)性。

4 結(jié)論

Lasso方法在變量選擇時具有很好的性質(zhì)，但是在處理海量的基因數(shù)據(jù)選擇相關(guān)致病基因時，會出現(xiàn)費時，消耗大的問題。于是為了更好地解決關(guān)聯(lián)基因位點的選擇，本文提出分而治之方法，此方法除了擁有Lasso的一般性質(zhì)外，還保留了強穩(wěn)定性，具有易排除偽相關(guān)變量的特性。實驗數(shù)據(jù)表明，分而治之方法通過了卡方與Fisher檢驗，并與Q-Q圖分析出的結(jié)果一致。但是，分而治之方法在K值的選擇上還需要進一步研究，以便得到最優(yōu)的K值。