基于遍歷基因組合的特征基因選取方法

李 杰

（大理學(xué)院數(shù)學(xué)與計(jì)算機(jī)學(xué)院，云南大理 671003）

隨著DNA基因芯片技術(shù)的提高，人們研究癌癥的方法也得到了快速發(fā)展。DNA 微陣列技術(shù)使得人們可以在基因水平上了解癌癥發(fā)生、癌變的機(jī)理。通常生物學(xué)上認(rèn)為癌癥是由于一個(gè)或者某幾個(gè)基因變異導(dǎo)致的?？茖W(xué)家們正致力于尋找產(chǎn)生癌變的基因或基因組合，近年來(lái)各種有效方法如雨后春筍般的涌現(xiàn)出來(lái)。

自1999年Golub〔1〕發(fā)表了關(guān)于急性白血病亞型基因芯片起，從基因芯片的角度研究各類疾病的人越來(lái)越多。由于基因芯片數(shù)據(jù)大部分是高維甚至超高維數(shù)據(jù)，而樣本量卻遠(yuǎn)遠(yuǎn)少于基因的維數(shù)。在統(tǒng)計(jì)學(xué)上，這樣的數(shù)據(jù)會(huì)導(dǎo)致“Curse of Dimensionality”和“Over-Fitting”問(wèn)題，統(tǒng)計(jì)學(xué)上解決這一問(wèn)題的方法通常有兩種：數(shù)據(jù)降維和變量選擇。數(shù)據(jù)降維方法是把高維矩陣在不損失信息或損失少量信息的情況下，把高維矩陣投影到一個(gè)低維矩陣的方法，進(jìn)而可以利用經(jīng)典的統(tǒng)計(jì)模型進(jìn)行分析，如 PCA〔2〕、SIR〔3〕、MAVE〔4〕等方法；變量選擇方法是通過(guò)對(duì)回歸系數(shù)的限制，在回歸過(guò)程中不斷把系數(shù)變?yōu)? 的變量剔除掉，直到達(dá)到所需要的變量個(gè)數(shù)為止，常用的方法有LASSO〔5〕、LAR〔6〕、ALASSO〔7〕、SCAD〔8〕、GLASSO〔9〕等。分類器可以用Logistic 回歸、支持向量機(jī)、決策樹C4.5和隨機(jī)森林的方法。

從可查的文獻(xiàn)顯示，選擇基因芯片中可能引發(fā)癌癥的基因，主要思路〔10〕如下：首先剔除無(wú)用和冗余基因，其次對(duì)剩下的基因進(jìn)行數(shù)據(jù)降維或者變量選擇，然后把選擇出的基因放到分類器中進(jìn)行分類，最后根據(jù)分類器的預(yù)測(cè)精度進(jìn)行修正和評(píng)價(jià)。但是第二步進(jìn)行變量選擇和數(shù)據(jù)降維時(shí)，并沒(méi)有考慮到基因之間的關(guān)系。很可能會(huì)出現(xiàn)這樣的情況：兩個(gè)基因單獨(dú)分析，對(duì)癌癥并沒(méi)有影響，但是合在一起對(duì)癌癥的產(chǎn)生卻有很大影響。就如雞蛋和糖精只能分開吃，不能一起吃的道理一樣。本文嘗試從這個(gè)角度出發(fā)，考慮存在此關(guān)系的基因?qū)Π┌Y的分類影響。

這方面的文獻(xiàn)不是很多。原因大概有3 點(diǎn)：①計(jì)算量大，基因個(gè)數(shù)通常以萬(wàn)計(jì)，不同的組合會(huì)呈幾何級(jí)數(shù)增長(zhǎng)；②計(jì)算受到了計(jì)算機(jī)內(nèi)存和計(jì)算速度的限制；③這種方法比較笨，效率低。但隨著計(jì)算機(jī)多核、多線程、分布式計(jì)算和云計(jì)算的快速發(fā)展，這類的問(wèn)題似乎可以嘗試運(yùn)行一下。

文章主要對(duì)結(jié)腸癌數(shù)據(jù)從考慮成對(duì)的基因的相互作用出發(fā)，用Logistic分類器對(duì)正常人和結(jié)腸癌患者進(jìn)行區(qū)分。結(jié)腸癌數(shù)據(jù)共包含62個(gè)樣本，其中22人為正常人的基因水平，40人為結(jié)腸癌的基因水平，共2 000 個(gè)基因?？紤]成對(duì)基因，則有199.9 萬(wàn)個(gè)不同的基因組合，遍歷所有可能的情況，根據(jù)預(yù)測(cè)精度和AIC準(zhǔn)則選擇最優(yōu)的基因組合并與已有的文獻(xiàn)進(jìn)行比較研究，最后根據(jù)前1 000個(gè)結(jié)果，畫出基因關(guān)系圖譜。

1 數(shù)據(jù)結(jié)構(gòu)及預(yù)處理

1.1 數(shù)據(jù)結(jié)構(gòu)結(jié)腸癌數(shù)據(jù)共包含62 個(gè)案例，其中22 個(gè)正常人，40 個(gè)癌癥患者；對(duì)每個(gè)患者運(yùn)用DNA 微陣列測(cè)得了2 000 個(gè)基因的表達(dá)水平。用Yi=1 表示第i個(gè)人是癌癥患者，Yi=0 表示第i個(gè)人是正常人，i=1,…,62，則Y是一個(gè)二元變量，每一個(gè)基因看成一個(gè)變量，分別用G1,G2,…,G2000表示。數(shù)學(xué)模型可以抽象為表1。

表1 結(jié)腸癌數(shù)據(jù)結(jié)構(gòu)

其中Xij，i=1,…,62,j=1,…,2 000 是第i個(gè)人第j個(gè)基因表達(dá)的對(duì)數(shù)值，該數(shù)據(jù)可以在BRB-Array-Tools主頁(yè)上下載。

1.2 預(yù)處理由于數(shù)據(jù)取對(duì)數(shù)后各個(gè)基因的表達(dá)水平還是有很大差距，因此需要按基因?qū)?shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，因此表1中的列需要做標(biāo)準(zhǔn)化變化，以便更好比較。一般在用分類器進(jìn)行分類之前，需要對(duì)冗余基因進(jìn)行處理，值得注意的是這2 000 個(gè)基因的表達(dá)水平都是非零，為了能準(zhǔn)確反映結(jié)果，沒(méi)有再剔除冗余基因。從2 000 個(gè)基因中任意選取2 個(gè)基因，不同的組合共有1 999 000 種，任意選擇3 個(gè)基因，不同的組合共有種。由于選擇3 個(gè)基因不同的情況達(dá)到了13.3 億，本文暫時(shí)不予考慮，而只考慮2個(gè)基因共同作用的情形。

2 模型方法

由于Y是一個(gè)分類的二元變量，在統(tǒng)計(jì)學(xué)中，二元離散變量的經(jīng)典模型是Logistic 模型，假設(shè)Y是因變量，x1,x2是自變量，Y取值為1 的概率為p，則Logistic模型為：

其中ε是隨機(jī)誤差項(xiàng)，β0，β1，β2是未知參數(shù)，可以通過(guò)迭代極大似然估計(jì)得到。由于只考慮2個(gè)基因聯(lián)合作用的情形，自變量的個(gè)數(shù)定位2 個(gè)。把上式進(jìn)行變換可以得到Y(jié)取1的概率為：

當(dāng)β0，β1，β2估計(jì)出來(lái)后，根據(jù)上式就可以計(jì)算第i個(gè)樣本取1的概率pi，當(dāng)pi＞0.5 時(shí)，把第i個(gè)樣本判定為癌癥患者，pi≤0.5 時(shí)判為正常人。根據(jù)預(yù)測(cè)值和實(shí)際值的比較，可以得到表2。

表2 預(yù)測(cè)值和實(shí)際值

其中n11表示實(shí)際為癌癥患者并且預(yù)測(cè)成癌癥患者的人數(shù)，n10表示實(shí)際為正常人但預(yù)測(cè)成癌癥患者的人數(shù)，n01表示實(shí)際為癌癥患者但預(yù)測(cè)成正常的人數(shù)，n00表示實(shí)際為正常人并且預(yù)測(cè)成正常的人數(shù)，n11+n10+n01+n00=n=62，預(yù)測(cè)精度可以表示為：

預(yù)測(cè)精度越大，效果越好。

3 評(píng)價(jià)標(biāo)準(zhǔn)

模型評(píng)價(jià)標(biāo)準(zhǔn)主要有兩部分組成：預(yù)測(cè)精度和統(tǒng)計(jì)指標(biāo)。第一部分主要看預(yù)測(cè)精度，預(yù)測(cè)精度越高越好，但是預(yù)測(cè)精度單個(gè)指標(biāo)并不能反映真實(shí)情況，甚至可能會(huì)出現(xiàn)錯(cuò)誤，例如有100個(gè)人，90個(gè)正常人，10 個(gè)病人，但是預(yù)測(cè)的結(jié)果把這100 個(gè)人全部預(yù)測(cè)成正常人，預(yù)測(cè)精度也達(dá)到了90%，但是它把病人全部預(yù)測(cè)成正常人，沒(méi)有一個(gè)預(yù)測(cè)正確，這種方法就沒(méi)有區(qū)分度，需要結(jié)合統(tǒng)計(jì)指標(biāo)。統(tǒng)計(jì)指標(biāo)本文采用了AIC準(zhǔn)則和β0，β1，β2的顯著性水平作為評(píng)價(jià)標(biāo)準(zhǔn)，在重要性排序上文獻(xiàn)中沒(méi)有統(tǒng)一說(shuō)法，根據(jù)需要擬定重要性大小按照預(yù)測(cè)精度，AIC信息和β0，β1，β2的顯著性水平依次遞減。

4 計(jì)算流程

計(jì)算流程共分為以下幾步：①標(biāo)準(zhǔn)化各個(gè)基因；②產(chǎn)生2個(gè)基因組合的所有可能結(jié)果；③遍歷所有可能的組合，在每一次遍歷過(guò)程中計(jì)算預(yù)測(cè)精度、AIC 和未知參數(shù)的檢驗(yàn)P值，并保存在txt 文檔中；④讀取最后的數(shù)據(jù)，按照重要性標(biāo)準(zhǔn)進(jìn)行排序，找出最合理的模型。流程圖見圖1。

圖1 計(jì)算流程圖

5 實(shí)驗(yàn)結(jié)果

本文采用的實(shí)驗(yàn)環(huán)境配置：Intel（R）Core（TM）i5-3470 CPU@3.20Hz，4GB 內(nèi)存，Windows 8.1 操作系統(tǒng)，Rstudio 平臺(tái)，用R 語(yǔ)言進(jìn)行編程。實(shí)驗(yàn)過(guò)程中對(duì)每個(gè)計(jì)算過(guò)程中的中間結(jié)果進(jìn)行保存，共得到了1 999 000行17列的數(shù)據(jù)文件，約225.3 MB。用R進(jìn)行分析得到的主要結(jié)果如下。

首先按照預(yù)測(cè)精確度對(duì)所有的數(shù)據(jù)進(jìn)行排序，在此基礎(chǔ)上再對(duì)AIC 從小到大進(jìn)行排序，同時(shí)檢驗(yàn)未知參數(shù)的t檢驗(yàn)P值是否小于0.05，結(jié)果顯示編號(hào)為X55187 和D14812 的基因組合得到最滿意的結(jié)果，預(yù)測(cè)精度為0.935 5，預(yù)測(cè)表格見表3。

表3 實(shí)驗(yàn)結(jié)果

AIC 為 29.671，在所有 199.9 萬(wàn)中 AIC 中排在第二位，最小的AIC為26.198。實(shí)際的模擬結(jié)果見表4。

表4 Logistic回歸結(jié)果

從表4中最后一列可以看出，包括截距項(xiàng)在內(nèi)，未知參數(shù)在置信水平0.01下都通過(guò)了檢驗(yàn)，參數(shù)值都是非零，具有統(tǒng)計(jì)學(xué)意義。根據(jù)表4，回歸方程可以寫為：

其中x1表示編號(hào)為X55187 的基因，x2表示編號(hào)為D14812 的基因。零假設(shè)下的離差為80.648 4，殘差為 23.671 1，廣義R2為：（80.6484-23.6711）/80.6484=0.7065。為了進(jìn)一步說(shuō)明該回歸結(jié)果的穩(wěn)定性，采用交叉留一檢驗(yàn)，即依次刪除掉62個(gè)樣本的1個(gè)樣本，剩下的61 樣本作為測(cè)試集，統(tǒng)計(jì)測(cè)試集的預(yù)測(cè)精度。見圖2。

圖2 交叉留一檢驗(yàn)預(yù)測(cè)精度趨勢(shì)圖

圖2顯示，預(yù)測(cè)精度大部分是0.934，其中有4次試驗(yàn)的預(yù)測(cè)精度達(dá)到了0.951，這說(shuō)明選擇出的基因組合具有較好的穩(wěn)健性。不會(huì)因?yàn)闃颖镜淖兓?，而?duì)預(yù)測(cè)精度產(chǎn)生較大波動(dòng)。綜上所述，編號(hào)為X55187和D14812的基因組合是比較適合區(qū)分結(jié)腸癌患者和正常人的。

因?yàn)楸闅v了所有可能情況，得到的結(jié)果豐富，為了進(jìn)一步展示其它結(jié)果，我們截取了199.9萬(wàn)中預(yù)測(cè)精度大于0.9的所有基因組合，共有640對(duì)基因組合，對(duì)這些基因組合進(jìn)行頻數(shù)分析，得到結(jié)果見表5。

表5 基因頻率分析結(jié)果

張靖等〔11〕研究結(jié)果顯示，基因R87126，H08393，M63391，X12671是引起結(jié)腸癌的可疑基因，跟本文中頻率分析結(jié)果相同的有兩個(gè)基因R87126，M63391，注意這里考察的是成對(duì)基因的頻率分析，也就是說(shuō)本文只是統(tǒng)計(jì)了各個(gè)基因出現(xiàn)的頻率，并沒(méi)有考慮成對(duì)基因之間內(nèi)部的關(guān)系。雖然R87126，M63391 在統(tǒng)計(jì)結(jié)果中出現(xiàn)的頻率最高，但是這兩者組合成的基因?qū)Σ](méi)有出現(xiàn)在640 對(duì)當(dāng)中，重新考察R87126，M63391組成的基因?qū)ΓY(jié)果顯示預(yù)測(cè)精度僅為87.10%，AIC為43.234，廣義R2為0.538 3，并且M63391基因的系數(shù)在0.05的置信水平下沒(méi)有通過(guò)檢驗(yàn)。因此并不是基因出現(xiàn)頻率越大，分類的預(yù)測(cè)精度就越高。

6 結(jié)論

本文在癌癥是由某個(gè)或者某幾個(gè)基因共同變異導(dǎo)致的假設(shè)下，采用遍歷所有基因組的方法，運(yùn)用Logistic 分類器，對(duì)結(jié)腸癌的基因進(jìn)行了篩選，按照預(yù)測(cè)精度和AIC 準(zhǔn)則得到了最優(yōu)的基因組合（X55187，D14812），并與已有文獻(xiàn)結(jié)果進(jìn)行比較，得出頻率高的基因組合預(yù)測(cè)精度并不一定高的結(jié)論。由于實(shí)驗(yàn)條件限制，本文中并沒(méi)有考察3 個(gè)以及3個(gè)以上基因組合對(duì)分類的影響。

〔1〕Alon U，Barkai N，Notterman D A，et al. Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays〔J〕.Proc Natl Acad Sci，1999，96：6745-6750.

〔2〕Pearson K. On lines and planes of closest fit to systems of points in space〔J〕. Philosophical Magazine，1901，2（6）：559-572.

〔3〕Li K C. Sliced inverse regression for dimension reduction〔J〕. Journal of the American Statistical Association，1991，86：316-327.

〔4〕Xia Yingcun，Tong Howell，Li W K.An adaptive estimation of dimension reduction space〔J〕. Journal of the Royal Statistical Society：Series B（Statistical Methodology），2002，64（3）：363-410.

〔5〕Tibshirani R. Regression shrinkage and selection via the lasso〔J〕.J Royal Statist Soc B，1996，58（1）：267-288.

〔6〕Zhou Hui. The Adaptive Lasso and Its Oracle Properties〔J〕. Journal of the American Statistical Association，2006，101（476）：1418-1429.

〔7〕Fan J，Li R. Variable selection via nonconcave penalized likelihood and its oracle properties〔J〕.J Amer Statist Assoc，2001，96：1348-1360.

〔8〕Noah S，Jerome F，Trevor H，et al.A sparse-group lass〔J〕.Journal of Computational and Graphical Statistics，2013，22（2）：231-245.

〔9〕Efron B，Hastie T J，Johnstone I，et al. Least Angle Regression〔J〕.Annals of Statistics，2004，32（2）：407-499.

〔10〕李杰，鄧麗君，唐勝男.基于BB—SIR方法的結(jié)腸癌特征基因提取〔J〕.世界科技研究與發(fā)展，2011，33（4）：588-591.

〔11〕張靖，胡學(xué)鋼，李培培，等.基于迭代Lasso 的腫瘤分類信息基因選擇方法研究〔J〕. 模式識(shí)別與人工智能，2014，27（1）：49-59.