張恒益，鄭惠玲

(西北農(nóng)林科技大學(xué) 動(dòng)物科技學(xué)院,陜西 楊陵712100)

生物醫(yī)學(xué)研究表明，許多疾病都有相應(yīng)的致病基因或易感基因。復(fù)雜疾病往往由基因變異引起，但致病基因的精確定位卻不易做到[1-2]。在組成DNA的堿基對(duì)中，一些特定位置的單核苷酸發(fā)生變異引起的DNA多態(tài)性稱為SNP (Single Nucleotide Polymorphism)位點(diǎn)。SNP是最常見(jiàn)的基因變異類型，識(shí)別致病SNP位點(diǎn)對(duì)精確定位致病基因具有重要意義。Pico等[3]通過(guò)識(shí)別致病SNP位點(diǎn)找到了孟德?tīng)柤膊〉娜毕莼?；Irina等[4]發(fā)現(xiàn)了克羅恩(Crohn)疾病相關(guān)的致病位點(diǎn)；Pernille等[5]發(fā)現(xiàn)了一些與情感信號(hào)傳導(dǎo)障礙相關(guān)的SNP位點(diǎn)。牛蜘蛛腿綜合征是一種以骨骼畸形為病理特征的先天致死性遺傳病，該病有兩個(gè)致病位點(diǎn)[6]。另外還有許多家畜遺傳疾病是由致病位點(diǎn)或基因引起，如肢蹄畸形、血友病、雞的卷羽、豬的應(yīng)激綜合征等[7]。

近年來(lái)，人們提出了許多識(shí)別致病SNP位點(diǎn)的方法[8-15]，其中性能優(yōu)異的邏輯斯蒂回歸及隨機(jī)森林算法目前被普遍使用[13]。邏輯斯蒂回歸是一種廣義線性回歸分析模型，解決了線性回歸模型不適用因變量為離散值的問(wèn)題，目前在疾病診斷、經(jīng)濟(jì)預(yù)測(cè)等領(lǐng)域應(yīng)用廣泛。隨機(jī)森林是一個(gè)包含多個(gè)決策樹的集成分類器，可顯著提升模型的精度及泛化性能。但是，以上兩種方法有一個(gè)共同的缺點(diǎn)是計(jì)算復(fù)雜度比較高，不適合實(shí)時(shí)處理數(shù)據(jù)量龐大的高維基因數(shù)據(jù)集。基于此，本文提出了利用K均值聚類算法識(shí)別遺傳疾病致病SNP位點(diǎn)的新算法，為實(shí)時(shí)處理大規(guī)模畜禽基因數(shù)據(jù)集提供參考。

1 材料與方法

1.1 數(shù)據(jù)

1.1.1 模擬數(shù)據(jù)集 利用文獻(xiàn)[16]提供的3種生成模擬數(shù)據(jù)集的模型。模型1具有位點(diǎn)乘法效應(yīng)；模型2具有規(guī)律的位點(diǎn)交互作用；模型3具有非規(guī)律的位點(diǎn)交互作用。因此，由模型3生成的模擬數(shù)據(jù)集相比模型1和模型2致病位點(diǎn)更難識(shí)別。為了更好地驗(yàn)證算法性能，本試驗(yàn)采用由模型3生成的模擬數(shù)據(jù)集。通過(guò)調(diào)節(jié)致病率參數(shù)，數(shù)據(jù)集中嵌入了兩個(gè)致病位點(diǎn)，編號(hào)為11和21。數(shù)據(jù)集由2 000個(gè)樣本組成，其中包括1 000個(gè)健康樣本和1 000個(gè)患病樣本，每個(gè)樣本包含1 000個(gè)位點(diǎn)。最小等位基因頻率(MAF)值取0.2，主效應(yīng)值取0.5，表示連鎖不平衡程度的值取0.33。模擬數(shù)據(jù)集的詳細(xì)情況如表1所示。

表1 模擬數(shù)據(jù)集詳細(xì)情況Table 1 The details of simulation data set

表2 真實(shí)數(shù)據(jù)集各位點(diǎn)數(shù)值編碼Table 2 The numerical coding of all SNPs of real data set

1.2 方法

1.2.1 邏輯斯蒂回歸算法 邏輯斯蒂回歸(Logistic Regression)實(shí)際上是一種分類模型，采用Sigmoid函數(shù)作為后驗(yàn)概率分布函數(shù)對(duì)輸入數(shù)據(jù)集進(jìn)行分類，常用于二分類問(wèn)題，由如下條件概率分布表示：

(1)

式中：X∈Rn表示輸入數(shù)據(jù)向量；W∈Rn表示系數(shù)權(quán)重向量；W·X表示二者的內(nèi)積；y∈{0,1}表示輸出類別。

1.2.2 隨機(jī)森林算法 2001年，Breiman和Cutler[18]基于決策樹(Decision Tree)提出了隨機(jī)森林 (Random Forests)算法。隨機(jī)森林本質(zhì)上是一種集成算法(Ensemble Learning)，屬于Bagging類型，即通過(guò)組合多個(gè)弱分類器(決策樹)的結(jié)果得到最終結(jié)果(投票),顯著提高了模型的精確度和泛化性能,原因就在于“隨機(jī)”和“森林”，一個(gè)使其具有抗過(guò)擬合能力，一個(gè)使其更加精準(zhǔn)。其工作原理如下：①?gòu)臄?shù)據(jù)集的總共m個(gè)特征中隨機(jī)選擇k(≤m)個(gè)特征，然后根據(jù)這k個(gè)特征建立決策樹；②重復(fù)選擇n次，建立n棵決策樹；③利用每棵決策樹進(jìn)行預(yù)測(cè)，得到n個(gè)預(yù)測(cè)結(jié)果；④計(jì)算每個(gè)預(yù)測(cè)結(jié)果的得票數(shù)，將得到票數(shù)最高的預(yù)測(cè)結(jié)果作為隨機(jī)森林算法的最終結(jié)果。

1.2.3 箱型圖統(tǒng)計(jì)法 箱型圖(box plot)是一種顯示一組數(shù)據(jù)分散情況的統(tǒng)計(jì)圖，能顯示出一組數(shù)據(jù)的上下限、中位數(shù)及上下四分位數(shù)，如圖1所示。

圖1 箱型圖Fig. 1 Box-plot

箱型圖可用于識(shí)別異常值，其識(shí)別溫和異常值及極端異常值的標(biāo)準(zhǔn)公式如下：

MA={x|xQ3+1.5I}

(2)

EA={x|xQ3+3I}

(3)

式中：Q1表示下四分位數(shù)；Q3表示上四位數(shù)；I=Q3-Q1；MA表示溫和異常值集合；EA表示極端異常值集合。相比普遍使用的基于正態(tài)分布的Pauta準(zhǔn)則，箱型圖不受極端值影響，不需要假定數(shù)據(jù)服從特定分布，異常值識(shí)別結(jié)果比較客觀[19]。為了保證致病SNP位點(diǎn)識(shí)別的準(zhǔn)確性，本文采用極端異常值標(biāo)準(zhǔn)公式(3)。

1.2.4 χ2檢驗(yàn) χ2檢驗(yàn)又稱卡方檢驗(yàn)[20]，其基本公式為：

(4)

式中：Ai為實(shí)際觀察頻數(shù)，n為樣本容量，Pi為理論概率，nPi為樣本理論頻數(shù)。

1.2.5 K均值聚類算法 K均值聚類算法(K-means clustering algorithm)是一種迭代求解的聚類分析算法，其步驟是隨機(jī)選取K個(gè)數(shù)據(jù)對(duì)象作為初始的聚類中心，然后計(jì)算每個(gè)對(duì)象與各個(gè)種子聚類中心之間的距離，把每個(gè)對(duì)象分配給距離其最近的聚類中心[21]。聚類中心以及分配給它們的對(duì)象代表一個(gè)聚類。每分配一個(gè)樣本，聚類的聚類中心會(huì)根據(jù)聚類中現(xiàn)有的對(duì)象被重新計(jì)算。這個(gè)過(guò)程不斷重復(fù)，不斷降低類簇的誤差平方和(Sum of Squared Error ，簡(jiǎn)記為SSE)，直到不再變化為止，得到最終結(jié)果。實(shí)際操作中終止條件通常取沒(méi)有(或最小數(shù)目)對(duì)象被重新分配給不同的聚類，或沒(méi)有(或最小數(shù)目)聚類中心再發(fā)生變化，或達(dá)到最大的迭代次數(shù)。

數(shù)據(jù)對(duì)象與聚類中心間的歐式距離計(jì)算公式為：

(5)

式中：X為數(shù)據(jù)對(duì)象；Ci為第 個(gè)聚類中心；m為數(shù)據(jù)對(duì)象的維度；xj和cij分別為X和Ci的第j個(gè)屬性。

整個(gè)數(shù)據(jù)集的誤差平方和SSE的計(jì)算公式為：

(6)

式中：k為類簇個(gè)數(shù)。

1.2.6 基于K均值聚類算法 基于K均值聚類算法、箱型圖統(tǒng)計(jì)法及卡方檢驗(yàn)，本文提出一種新的識(shí)別致病SNP位點(diǎn)的算法。該算法工作原理如下：①利用K均值聚類算法將每個(gè)位點(diǎn)的所有樣本聚為2類，即健康樣本和患病樣本。②根據(jù)數(shù)據(jù)文件提供的樣本類別信息，按照公式(7)計(jì)算聚類正確率。聚類正確率越高的位點(diǎn)說(shuō)明區(qū)分健康樣本和患病樣本的能力越強(qiáng)，是致病位點(diǎn)的可能性越大；反之，是致病位點(diǎn)的可能性越小。

(7)

式中：N表示樣本總數(shù)；TP表示聚類正確的健康樣本；TN表示聚類正確的患病樣本。

③利用箱型圖統(tǒng)計(jì)法篩選異常值，即致病SNP位點(diǎn)。

④利用χ2檢驗(yàn)進(jìn)一步校驗(yàn)箱型圖篩選出的致病SNP位點(diǎn)。

1.2.7 程序和運(yùn)行環(huán)境 為了驗(yàn)證所提出的基于K均值聚類算法的性能，將該算法與當(dāng)前廣泛應(yīng)用的邏輯斯蒂回歸和隨機(jī)森林兩個(gè)主流算法在兩個(gè)數(shù)據(jù)集上進(jìn)行識(shí)別準(zhǔn)確性及運(yùn)行時(shí)間的比較。三種算法程序均采用MATLAB語(yǔ)言編寫，其中本文編寫了基于K均值聚類算法及邏輯斯蒂回歸算法，分別調(diào)用了實(shí)現(xiàn)K均值聚類的kmeans和實(shí)現(xiàn)邏輯斯蒂回歸的glmfit兩個(gè)MATLAB內(nèi)置函數(shù)；隨機(jī)森林算法程序來(lái)自GitHub MATLAB程序包，參數(shù)值取缺省值。所有試驗(yàn)均使用64位Win10系統(tǒng)，在MATLAB 2016軟件上進(jìn)行，其中CPU為Intel(R)Core (TM) i7-4790 ( 3.6 GHz)，內(nèi)存大小為8GB。

2 結(jié)果與分析

2.1 基于K均值聚類算法

2.1.1 模擬數(shù)據(jù)集 ① 聚類結(jié)果。利用K均值聚類算法對(duì)模擬數(shù)據(jù)集每個(gè)位點(diǎn)進(jìn)行聚類并計(jì)算正確率。正確率排名前十的位點(diǎn)及其卡方檢驗(yàn)統(tǒng)計(jì)量值如表3所示。各個(gè)位點(diǎn)的聚類正確率如圖2所示。

表3 模擬數(shù)據(jù)集聚類正確率排名前10位點(diǎn)的計(jì)算結(jié)果Table 3 The details of top 10 SNPs of clustering accuracyin simulation data set

由表3可得，位點(diǎn)SNP11的聚類正確率比SNP21高0.09%，而比位點(diǎn)SNP185高3.43%，因此可以初步推斷，模擬數(shù)據(jù)集位點(diǎn)SNP11和SNP21為疑似致病位點(diǎn)。

②致病位點(diǎn)篩選。利用箱型圖統(tǒng)計(jì)法公式(3)篩選聚類正確率的極端異常值，即致病位點(diǎn)。由各個(gè)位點(diǎn)的聚類正確率值可以計(jì)算出箱型圖極端異常值閾值上限為0.5385。由圖2和表3可以篩選出致病位點(diǎn)為SNP11和SNP21。

圖2 模擬數(shù)據(jù)集每個(gè)位點(diǎn)聚類正確率Fig. 2 Clustering accuracy of eachSNP in simulation data set

③卡方檢驗(yàn)驗(yàn)證。利用卡方檢驗(yàn)對(duì)箱型圖篩選結(jié)果做進(jìn)一步的驗(yàn)證。由公式(4)可以計(jì)算出每個(gè)位點(diǎn)的卡方值，如圖3所示。

由于數(shù)據(jù)分為2類和3種基因型，所以卡方檢驗(yàn)自由度為(2-1)×(3-1)=2。

Bonferroni原理[22]：如果在同一數(shù)據(jù)集上同時(shí)檢驗(yàn)n個(gè)獨(dú)立的假設(shè)，那么用于每一假設(shè)的統(tǒng)計(jì)顯著水平應(yīng)為僅檢驗(yàn)一個(gè)假設(shè)時(shí)顯著水平的1/n。

由于模擬數(shù)據(jù)集含有1 000個(gè)位點(diǎn)，所以顯著水平取值為0.05/1000=5×10-5,查表可得卡方檢驗(yàn)臨界值為19.80。

圖3 模擬數(shù)據(jù)集每個(gè)位點(diǎn)卡方值Fig. 3 Chi square test statistic value ofeach SNP in simulation data set

由表3和圖3可知，只有位點(diǎn)SNP11和SNP21的卡方值高于臨界值19.80，這與本文提出的基于K均值聚類算法篩選結(jié)果一致。

2.1.2 真實(shí)數(shù)據(jù)集 ①聚類結(jié)果。利用K均值聚類算法對(duì)真實(shí)數(shù)據(jù)集每個(gè)位點(diǎn)進(jìn)行聚類然后計(jì)算聚類正確率。正確率排名前十的位點(diǎn)及其卡方檢驗(yàn)統(tǒng)計(jì)量值見(jiàn)表4，所有位點(diǎn)的聚類正確率見(jiàn)圖4。

表4 聚類正確率排名前10位點(diǎn)的計(jì)算結(jié)果Table 4 The details of top 10 SNPs of clusteringaccuracy of real data set

圖4 真實(shí)數(shù)據(jù)集每個(gè)位點(diǎn)聚類正確率Fig. 4 Clustering accuracy of each SNP in real data set

從表4可得，真實(shí)數(shù)據(jù)集SNP2938位點(diǎn)的聚類正確率比SNP962高4.46%，而SNP962比SNP1541高0.18%，因此可以初步推斷，真實(shí)數(shù)據(jù)集中位點(diǎn)SNP2938是疑似致病位點(diǎn)，其相應(yīng)的位點(diǎn)名稱為rs2273298。

②致病位點(diǎn)篩選。利用箱型圖統(tǒng)計(jì)法公式(3)篩選聚類正確率的極端異常值。由各個(gè)位點(diǎn)的聚類正確率值可以計(jì)算出箱型圖極端異常值閾值上限為0.5590。由圖4和表4可以篩選致病位點(diǎn)為SNP2938，即rs2273298。

③卡方檢驗(yàn)驗(yàn)證。利用卡方檢驗(yàn)對(duì)箱型圖篩選結(jié)果做進(jìn)一步的驗(yàn)證。由公式(4)可以計(jì)算出每個(gè)位點(diǎn)的卡方值，如圖5所示。由于真實(shí)數(shù)據(jù)集依然分為2類和3種基因型，所以卡方檢驗(yàn)自由度依然為(2-1)×(3-1)=2。

依據(jù)Bonferroni原理，由于真實(shí)數(shù)據(jù)集含有9 445個(gè)位點(diǎn)，所以顯著水平取值為0.05/9445=5.2938×10-6，查表可得卡方檢驗(yàn)臨界值為24.2979。

圖5 真實(shí)數(shù)據(jù)集每個(gè)位點(diǎn)卡方值Fig. 5 Chi-square test statistic value of eachSNP in real data set

由表4和圖5很容易看出，只有位點(diǎn)SNP2938的卡方值高于臨界值，這與本文提出的基于K均值聚類算法篩選結(jié)果一致。

2.2 邏輯斯蒂回歸算法和隨機(jī)森林算法

2.2.1 模擬數(shù)據(jù)集 邏輯斯蒂回歸算法和隨機(jī)森林算法通常被用于致病位點(diǎn)識(shí)別[14,22-23],兩種算法都選擇分類正確率作為篩選致病位點(diǎn)的指標(biāo)。分類正確率的計(jì)算方法類似公式(7)，只需將聚類換成分類即可。同樣使用箱型圖統(tǒng)計(jì)法篩選致病位點(diǎn)。經(jīng)計(jì)算，邏輯斯蒂回歸算法的極端異常值閾值上限為0.5375，隨機(jī)森林算法的極端異常值閾值上限為0.5395。其中，邏輯斯蒂回歸分析函數(shù)glmfit返回參數(shù)中包括P值。同樣，依據(jù)Bonferroni原理，模擬數(shù)據(jù)集的顯著性水平依然取值為0.05/1000=5×10-5。邏輯斯蒂回歸算法分類正確率排名前十的位點(diǎn)及相應(yīng)P值見(jiàn)表5。

表5 邏輯斯蒂回歸算法分類準(zhǔn)確率排名前10位點(diǎn)及P值Table 5 The top 10 SNPs using Logistic Regression

由表5可以看出，邏輯斯蒂回歸算法準(zhǔn)確識(shí)別出了兩個(gè)致病位點(diǎn)，這兩個(gè)致病位點(diǎn)也都通過(guò)了P值檢驗(yàn)。

隨機(jī)森林算法分類正確率排名前十的位點(diǎn)見(jiàn)表6。由表6可以看出，隨機(jī)森林算法也準(zhǔn)確識(shí)別出了模擬數(shù)據(jù)集的致病位點(diǎn)。

表6 隨機(jī)森林算法分類正確率排名前10位點(diǎn)Table 6 The top 10 SNPs of classification accuracyusing Random Forest

2.2.2 真實(shí)數(shù)據(jù)集 方法和步驟與模擬數(shù)據(jù)集一樣。經(jīng)計(jì)算，邏輯斯蒂回歸算法的極端異常值閾值上限為0.5590，隨機(jī)森林算法的極端異常值閾值上限為0.5600。同樣，依據(jù)Bonferroni原理，真實(shí)數(shù)據(jù)集的顯著性水平依然取值為0.05/9445=5.2938×10-6。邏輯斯蒂回歸算法分類正確率排名前十的位點(diǎn)及相應(yīng)P值見(jiàn)表7。由表7可以看出，邏輯斯蒂回歸對(duì)真實(shí)數(shù)據(jù)集也正確識(shí)別出了致病位點(diǎn)，這個(gè)致病位點(diǎn)也通過(guò)了P值檢驗(yàn)。

表7 邏輯斯蒂回歸算法分類準(zhǔn)確率排名前10位點(diǎn)及P值Table 7 The top 10 SNPs using Logistic Regression

隨機(jī)森林算法分類正確率排名前十的位點(diǎn)見(jiàn)表8。由表8可以看出，隨機(jī)森林算法也準(zhǔn)確識(shí)別出了致病位點(diǎn)。

表8 隨機(jī)森林算法分類正確率排名前10位點(diǎn)Table 8 The top 10 SNPs of classification accuracyusing Random Forest

2.3 三種算法程序運(yùn)行時(shí)間比較

三種算法的運(yùn)行時(shí)間如表9所示，其中本文提出的基于K均值聚類算法和邏輯斯蒂回歸算法的數(shù)據(jù)都是運(yùn)行10次取平均值的結(jié)果，而隨機(jī)森林算法數(shù)據(jù)是運(yùn)行一次的結(jié)果，因?yàn)殡S機(jī)森林算法運(yùn)行時(shí)間太長(zhǎng)。

表9 三種算法運(yùn)行時(shí)間對(duì)比Table 9 Running time comparison of three algorithms s

由表9可以看出，本文提出的基于K均值聚類算法用于模擬數(shù)據(jù)集和真實(shí)數(shù)據(jù)集的運(yùn)行時(shí)間相比其它兩種算法都有很大的優(yōu)勢(shì)。

3 討 論

本研究在對(duì)遺傳位點(diǎn)進(jìn)行數(shù)值編碼的基礎(chǔ)上，使用K均值聚類算法，結(jié)合箱型圖與卡方檢驗(yàn)，提出了基于K均值聚類算法的一種識(shí)別致病SNP位點(diǎn)的新方法。該方法首先利用K均值聚類計(jì)算各位點(diǎn)的聚類正確率，然后應(yīng)用箱型圖統(tǒng)計(jì)法篩選出致病SNP位點(diǎn)，再用Bonferroni原理和卡方檢驗(yàn)對(duì)結(jié)果做進(jìn)一步的校驗(yàn)。最后與邏輯斯蒂回歸算法和隨機(jī)森林算法進(jìn)行了比較。

在識(shí)別結(jié)果方面，三種算法均準(zhǔn)確識(shí)別出了致病位點(diǎn)，與方雅蘭[14]、陳劍[22]和王宇琛等[23]研究結(jié)果一致。在閾值方面，對(duì)模擬數(shù)據(jù)集，隨機(jī)森林算法閾值優(yōu)于本文提出的基于K均值聚類算法及邏輯斯蒂回歸算法；對(duì)真實(shí)數(shù)據(jù)集，本文提出的基于K均值聚類算法及邏輯斯蒂回歸算法閾值優(yōu)于隨機(jī)森林。在運(yùn)行時(shí)間方面，對(duì)模擬數(shù)據(jù)集，邏輯斯蒂回歸算法大約是本文提出的基于K均值聚類算法4倍；對(duì)真實(shí)數(shù)據(jù)集，邏輯斯蒂回歸算法大約是本文提出的基于K均值聚類算法12倍。對(duì)模擬數(shù)據(jù)集，隨機(jī)森林算法大約是本文提出的基于K均值聚類算法670倍；對(duì)真實(shí)數(shù)據(jù)集，隨機(jī)森林算法大約是本文提出的基于K均值聚類算法1000倍。因此，在運(yùn)行速度方面，相比邏輯斯蒂回歸算法和隨機(jī)森林算法，本文提出的基于K均值聚類算法優(yōu)勢(shì)巨大，并且數(shù)據(jù)集規(guī)模越大優(yōu)勢(shì)越大。K均值聚類算法優(yōu)勢(shì)巨大的原因是：K均值聚類算法的復(fù)雜度為O(N)，而邏輯斯蒂回歸算法的復(fù)雜度為O(NlogN)，隨機(jī)森林算法的復(fù)雜度為O(MNlogN)，其中N為樣本數(shù)，M為樹的棵數(shù)。

本文提出的基于K均值聚類算法不但可以識(shí)別致病位點(diǎn)，也可以識(shí)別致病基因，借助其速度優(yōu)勢(shì)可用于畜禽遺傳疾病的大規(guī)模檢測(cè)，對(duì)預(yù)防遺傳疾病在畜禽中的散播，凈化譜系，提高選種準(zhǔn)確性，加快畜禽育種進(jìn)展，進(jìn)而提高畜牧業(yè)生產(chǎn)水平具有重要意義。