劉月麗，覃錫忠，賀三剛，李文蓉，王 悅，賈振紅，劉明軍

LIU Yueli1,QIN Xizhong1,HE Sangang2,LI Wenrong2,WANG Yue1,JIAZhenhong1,LIU Mingjun2

1.新疆大學(xué) 信息科學(xué)與工程學(xué)院，烏魯木齊 830046

2.新疆畜牧科學(xué)院 生物技術(shù)研究所，農(nóng)業(yè)部草食家畜繁育生物技術(shù)重點開放實驗室，新疆維吾爾自治區(qū)動物生物技術(shù)重點實驗室，烏魯木齊 830046

1.College of Information Science and Engineering,Xinjiang University,Urumqi 830046,China

2.The Key Laboratory of Livestock Reproduction&Breed Biotechnology of MOA,The Key Laboratory of Animal Biotechnology of Xinjiang,Xinjiang Academy of Biotechnological Center,Urumqi 830046,China

1 引言

生物品種鑒別有著廣泛和重要的應(yīng)用價值。品種鑒別促進了對遺傳信息的有效管理[1]；為育種策略的制定與實施奠定良好基礎(chǔ)[2]；為生物品牌產(chǎn)品的認證提供了有效信息[3]；更為解決食品安全問題開辟了新途徑[4]。在以往品種鑒別研究中少有涉及羊的種類鑒別，而我國有著飼養(yǎng)羊群的悠久歷史，擁有豐富的羊群遺傳資源，據(jù)《中國畜禽遺傳資源志·羊志》記載的我國地方羊品種就達42個[5]，約占世界羊品種的5%[1]。本文將SNP位點應(yīng)用于羊的品種分類。

在早期的品種鑒別研究中，利用微衛(wèi)星和ALFP標記進行動物鑒別[6-7]。隨著基因芯片技術(shù)的飛速發(fā)展，SNP位點漸漸用于動物品種及其副產(chǎn)品的鑒別中。與微衛(wèi)星標記相比SNP位點的優(yōu)勢在于：基因分型錯誤率低，在基因組中分布廣泛，檢測速度快，易于標準化[8]。品種鑒別是對基因數(shù)據(jù)實際應(yīng)用的一次探索性研究。已有的相關(guān)研究中將遺傳信息與基本的數(shù)學(xué)統(tǒng)計方法相結(jié)合用于品種鑒別：Pfaff等人利用δ方法[9]，以兩個物種間的等位基因頻率絕對差為判別標準進行分類；Weir等人利用 Wright’s FST方法[10]，依賴于預(yù)先定義的兩個物種間的等位基因頻率的差異最大化進行判別。然而δ和Wright’s FST只可以用于兩個種群的判別，并且沒有清晰的統(tǒng)計特性定義。為解決兩個品種以上的判別，Rosenberg等人[11]提出了一種相關(guān)性衡量的方法，使用互信息（In）描述相關(guān)性，以此來表示不同品種的FST之間的關(guān)系。很明顯，上述的傳統(tǒng)方法都沒有考慮到SNP數(shù)據(jù)的高維小樣本的特點[12]（本文中維度就達到了近50000），且都需要計算等位基因頻率，大大增加了計算復(fù)雜度和運算量。為此，Paschou等人[13]提出利用PCA提取SNP位點的方法。

PCA雖然可以有效地降低樣本維數(shù)，但是對于品種鑒別來說SNP位點個數(shù)還是過多。而且PCA對非線性數(shù)據(jù)的分類效果并不好，也無法評估所選SNP位點的重要程度。為了有效解決這些問題，本文將PCA和隨機森林相結(jié)合提取高信息量的SNP位點。隨機森林有以下優(yōu)點，可以利用數(shù)據(jù)之間相互依賴的信息，構(gòu)建多元分類規(guī)則；在每個類別邊界上有很好的非線性；提供了變量重要性的衡量方法。目前還未被用于品種鑒別的高信息量SNP位點的識別中。

為驗證方法的可行性，本文針對6種綿羊llumina OvineSNP50的SNP數(shù)（AustralianPollDorset羊、AustralianSuffolk羊、MilkLacaune羊、Soay羊、德國肉用美利奴羊GMS、哈薩克羊KSK），先利用PCA進行降維，去除冗余SNP位點，得到一個特征SNP子集。再利用隨機森林評估這些SNP位點的重要性，并選出重要性排名靠前的SNP位點用于分類。最后，通過分類測試實驗，驗證了本文方法能夠利用少量高信息量的SNP將未知品種的羊正確分類。

2 PCA和隨機森林相結(jié)合提取SNP位點

PCA和隨機森林相結(jié)合篩選高信息量SNP位點算法流程如圖1。

2.1 SNP數(shù)據(jù)集與預(yù)處理

先將SNP基因型數(shù)據(jù)進行編碼，但是SNP位點采樣時不可避免地存在缺失值，因此對整個數(shù)據(jù)集進行缺失值插補。插補概率遵循哈溫伯格平衡定律，使得SNP位點在不同品種中分布更加均衡，不會引入人為的偏倚。最后將數(shù)據(jù)分為訓(xùn)練集和測試集。

2.2 基于PCA的SNP位點提取

PCA是一種線性降維技術(shù)，可以從過于“豐富”的數(shù)據(jù)信息中捕獲最重要的元素和結(jié)構(gòu)。在遺傳學(xué)中每個個體可以提取出來數(shù)以萬記的SNP位點。由于維數(shù)過高，需要先將這些數(shù)據(jù)進行編碼，再利用PCA和奇異值分解（Singular Value Decomposition，SVD）提取主SNP位點。

將編碼后的SNP數(shù)據(jù)看成矩陣Xm×n，其中m為樣本個數(shù)，n為SNP位點數(shù)目；再用SVD對矩陣X進行運算，返回m個正交向量ui,n個正交向量vi以及m個非負的奇異值δi。矩陣X表示成如下形式：

其中，左奇異向量ui與X的列（SNP位點）有關(guān)，是X矩陣列向量的線性組合，稱作特征SNP[14]?？梢员Ａ羯贁?shù)特征SNP，達到降維的目的。再利用所選擇的特征SNP表示所有樣本，進行下一步的數(shù)據(jù)分析（如聚類或者分類）。

然而特征SNP是一種數(shù)學(xué)抽象，并不代表實際的SNP位點。Drineas等人[14-16]證明，由特征SNP張成的子空間所得到的奇異向量，可以近似表示矩陣X中的奇異向量?；谶@一理論，所選出來的特征SNP就和實際SNP位點有關(guān)。假設(shè)有k個主成分，因此有k個特征SNP，則原始矩陣的列可以用前k個特征SNP張成的子空間中的列近似表示，如下式：

式中，表示第j個元素的第i個右奇異向量。X第j列由所有左奇異向量和其相應(yīng)的奇異值組成，其中δiui是這個線性組合的系數(shù)。

圖1 PCA與隨機森林相結(jié)合的算法流程圖

最后，利用對X矩陣的列進行評分，見式（3），選出最大的列組成特征SNP子集。

Drineas等人[15-16]證明，從X矩陣中，每次以概率pj獨立隨機選取j列，則被選擇的這些列的前k個左奇異向量和原始矩陣的前k個左奇異向量十分接近。在實際應(yīng)用中，選擇pj值最大的列就可以達到很好的效果。

2.3 用于分類的隨機森林算法

在實際應(yīng)用中經(jīng)過PCA預(yù)選后的SNP位點數(shù)目仍然很多，無法評估所選SNP位點的重要程度。因此需要建立一個規(guī)則，將樣本分配到對應(yīng)的品種中。出于這些目的本文利用隨機森林進行分類。

隨機森林是一種集成分類器[17]，在建模時有兩個重要參數(shù)。一個是隨機森林中決策樹的數(shù)目，另一個是樹節(jié)點預(yù)選變量的個數(shù)。第一個參數(shù)決定了整片隨機森林的規(guī)模，第二個參數(shù)決定了單棵決策樹的情況。算法流程如圖2所示。圖中每個節(jié)點都基于單一特征進行分類，每個分支的結(jié)束為終端節(jié)點。終端節(jié)點根據(jù)決策樹的路徑對樣本的類別進行預(yù)測。終端節(jié)點的顏色表示預(yù)測的類別。對于樣本的最終預(yù)測結(jié)果是基于所有的決策樹的預(yù)測結(jié)果。

圖2 隨機森林的流程圖（顯示兩棵樹的細節(jié)）

隨機森林中每棵決策樹的訓(xùn)練集都是從原始訓(xùn)練集中抽取m個樣本，每個樣本的容量和原始訓(xùn)練集是一樣的。這種方法叫作自助抽樣法（Bootstrap Sampling）。這樣每次有30%～40%的數(shù)據(jù)沒被抽到，這些數(shù)據(jù)對于每棵樹來說就是“Out-of-Bag（OOB）”數(shù)據(jù)。在隨機森林的生成過程中利用這些數(shù)據(jù)在內(nèi)部進行評估，可以得到相應(yīng)的OOB錯誤率[18]，因此本文利用OOB錯誤率來評估隨機森林的分類效果。如果OOB錯誤率為0，就表明每個樣本都被正確分類了。

2.4 位點重要性評價準則

隨機森林中每棵樹的節(jié)點上，隨機抽取d個特征用于二分類，并且d＜＜D,D是原始數(shù)據(jù)集中特征值的個數(shù)。父節(jié)點np可以根據(jù)Gini指數(shù)來劃分子節(jié)點nl和nr。節(jié)點n的Gini指數(shù)可以用式（4）表示：

其中，pc代表在節(jié)點n處樣本屬于c類的相對概率。為實現(xiàn)最佳的二分類，兩個子節(jié)點之間的Gini指數(shù)差值要取最大，Gini指數(shù)差值計算見式（5）：

其中，pl和pr分別表示np分裂為子節(jié)點nl和nr的比例值。

Gini指數(shù)可以評價的重要程度。還有一種衡量標準——平均精度下降（Mean Accuracy Decrease），可以作為節(jié)點分裂的標準，同時也可以評價用于分類的變量的重要程度。

3 實驗及結(jié)果分析

3.1 SNP數(shù)據(jù)集與預(yù)處理

本文SNP位點數(shù)據(jù)來自于6個品種的羊，共計696個樣本，其中Australian Poll Dorset 108個，Australian Suffolk 109個，Milk Lacaune 103個，Soay 109個，GMS 106個，KSK 161 個。Australian Poll Dorset、Australian Suffolk、Soay、Milk Lacaune來自于公開數(shù)據(jù)源ISGC（International Sheep Genomic Consortium），每個樣本SNP位點數(shù)目46013個；而GMS和KSK由新疆畜牧科學(xué)院提供。共計26條常染色體，每個樣本SNP位點數(shù)目54241個，包含了ISGC的所有位點；兩類樣本合并處理，用于實驗的SNP位點個數(shù)共計46013個。然后進行數(shù)據(jù)預(yù)處理，過程如下：

（1）SNP位點編碼?；蛐虯A編碼為0；基因型AB編碼為1；基因型BB編碼為2，編碼后的SNP位點數(shù)據(jù)構(gòu)成矩陣Xm×n,m=696,n=46013。

（2）缺失值插補。利用0、1、2對整個數(shù)據(jù)集進行隨即插補。

（3）訓(xùn)練集與測試集。從每個品種中隨機抽取30%的樣本組成測試集（Australian Poll Dorset 32只，Australian Suffolk 32只，MilkLacaune 30只，Soay 32只，GMS 31只，KSK 48只），用于測試最后的品種分配。訓(xùn)練集由余下的70%樣本數(shù)據(jù)組成。

3.2 利用PCA選擇高信息量的SNP子集

本文使用R語言3.3.0版本進行PCA計算。為了減少計算時間，分別按照每一條染色體進行計算，提高計算效率。每一個主成分所提取的信息量用方差來度量，其中某一個主成分方差的貢獻就等于原指標相關(guān)矩陣相應(yīng)的特征值δi，則第i個主成分的方差貢獻率見式（6）：

Ti值越大，說明相應(yīng)的主成分反映綜合信息的能力越強。計算后的PCA方差貢獻率見圖3。

圖3 所有SNP位點在前10個主成分上的方差貢獻率

為保證準分類的準確性，最重要的是所選主成分是否涵蓋了不同品種的差異，而不是累計方差貢獻率的百分比。如圖3所示，在PC2和PC3之間，主成分的方差貢獻率值有大幅度的下降，因此保留前兩個主成分，用于減少SNP位點的數(shù)目。再根據(jù)式（3）計算每個SNP標記的得分。利用得分對每條染色體上的SNP位點進行排名，選擇每條染色體上排名前20的SNP位點，組成含有520個SNP位點的縮減集。

經(jīng)過PCA計算后，SNP位點個數(shù)從46013減少到了520個。圖4、圖5展示了從520個SNP位點的協(xié)方差矩陣中分別提取出前兩個和前3個主成分構(gòu)成的空間。對于保留前3個主成分僅僅是以探索為目的，和Paschou等人[13]提出的基于PCA的甄選算法不同，本文對PCA結(jié)果繪圖僅用于評價減少SNP的數(shù)量是否會丟失能夠表達品種差異的相關(guān)信息。

圖4 520個SNP位點的前兩個主成分分析圖

圖5 520個SNP位點的前3個主成分分析圖

由圖4可以看出，在由前兩個主成分組成的空間中GMS和KSK以及MilkLacaune之間有稍微的混疊，但是在三維空間圖5中GMS和KSK以及MilkLacaune之間的混疊就消失了。PCA預(yù)先選擇出的520個位點可以很好地表示不同品種之間的差異。

3.3 隨機森林重要性評價及分類模型建立

利用隨機森林，以預(yù)選出的520個SNP位點為基礎(chǔ)，根據(jù)兩種排名方式（Gini指數(shù)下降和平均精度下降），選擇具有更大差異的SNP位點。根據(jù)實際應(yīng)用保留排名前48和排名前96的SNP位點，用于最后的預(yù)測。對于這4種要求（兩種排名方式×兩種SNP集合），訓(xùn)練相應(yīng)的隨機森林模型，并計算對應(yīng)的OOB錯誤率。最后，用測試集來測試最終的分類效果。

上節(jié)提到構(gòu)建隨機森林模型有兩個重要參數(shù)，繪制圖6（a）隨機森林中決策樹的個數(shù)與相關(guān)誤差的關(guān)系，圖6（b）樹節(jié)點預(yù)選的變量個數(shù)與相關(guān)誤差的關(guān)系，通過這兩張圖來選擇參數(shù)值。

圖6 隨機森林參數(shù)與錯誤率之間的關(guān)系

由圖6（a）可知：該模型的決策樹數(shù)量在200以內(nèi)時，模型誤差會出現(xiàn)較大的波動。當(dāng)決策樹數(shù)量大于200后，模型逐漸趨于穩(wěn)定，但還是有少許波動。通過觀察發(fā)現(xiàn)在決策樹數(shù)量為600時該模型的錯誤率最低，因此選擇500為隨機森林中決策樹的數(shù)量。

由圖6（b）可知：該模型在5之前錯誤率較高，在Mtry=6處錯誤率最低，之后錯誤率又有升高，因此每棵樹節(jié)點預(yù)選變量的個數(shù)設(shè)置為6。

分別按照Gini指數(shù)下降和平均精度下降選出的SNP位點有一定的重疊。排名前48的SNP中有35個相同，排名前96的里面有68個相同。在Boulesteix等人[19]的研究中提出，基于Gini指數(shù)的排名靠前的位點的MAF（最小等位基因頻率）值都較大。這也體現(xiàn)了本文分析方法的特性：剔除了基因頻率較小的SNP位點。這樣，實驗中篩選出的SNP位點就能夠捕捉到每個品種的普遍特性，也就是說所選出的SNP位點對于這幾個品種來說具有較好的普適性。48個SNP位點的Gini指數(shù)如圖7所示，可以明顯看出不同的SNP位點對于分類的貢獻不同，更利于對貢獻較大的SNP位點進一步研究。

圖7 48個SNP位點的Gini指數(shù)值

基于Gini指數(shù)選出的兩種SNP集合的位點分布情況如圖8所示。由圖可知，在10號染色體上兩種SNP集合（48和96）都有較多的SNP位點被選中（48中有6個，96中也有6個）。在96個SNP中，被選中位點最多的是1號染色體，有8個SNP位點。在48個SNP中有7個染色體上沒有被選中的位點，分別是4、7、8、9、14、18、21號染色體。并且在96面板中，任何一條染色體上被選中的SNP位點個數(shù)都沒有達到20個。說明本文選擇出了每條染色體上最具代表性的SNP位點，也說明PCA預(yù)選出來的每條染色體上排名前20的SNP位點，是可以捕獲到高信息量的SNP位點（用于品種分類）。

圖8 被選中的SNP位點在每個染色體上的分布情況

3.4 隨機森林的分類結(jié)果評價

為驗證PCA可預(yù)選SNP位點，隨機森林可以進一步減少用于分類的SNP位點個數(shù)并達到較好的準確率。基于4種數(shù)據(jù)集分別隨機抽取48和96個SNP位點對隨機森林的分類結(jié)果進行測試，結(jié)果見表1。

表1 利用隨機森林的分類結(jié)果

方式1：從46013個SNP位點中隨機抽取48個和96個SNP位點，然后用訓(xùn)練好的隨機森林模型進行分類。OOB錯誤率達到了12.02%，而隨機選取96個SNP位點在進行分類時，降到了8.96%，說明保留的位點數(shù)越多其準確率也就越高。

方式2：從PCA預(yù)選出的520個SNP位點中隨機抽取48個和96個SNP位點，再利用隨機森林分類?？梢钥吹絆OB錯誤率有明顯的降低，這也進一步說明，利用PCA評分選出的SNP位點與隨機選取的SNP位點相比更具價值。

方式3：根據(jù)Gini指數(shù)減少的排名選取48和96個SNP位點。方式4：根據(jù)平均精度減少的排名選取48和96個SNP位點。OOB錯誤率相對于僅僅根據(jù)PCA評分選出的SNP位點來說，利用隨機森林兩種排名方式選出的SNP位點更具有價值能夠進一步降低錯誤率，提高計算效率。

不同的分類方法與本文方法對比結(jié)果見表2。

可以看到相對于文獻[20-21]，本文中每個品種的參考樣本是較多的，這表明本文方法能夠捕獲到更多的品種差異信息，會降低由冗余SNP位點導(dǎo)致錯分的可能性[22]。在文獻[20]中不同品種的分類需要根據(jù)決策樹的不同層次建立不同大小的SNP集合，這樣才能夠保證到達要求的準確性，這個過程較為復(fù)雜。本文的隨機森林利用相似結(jié)構(gòu)的決策樹，解決了不同層次的分類問題，優(yōu)化計算過程，節(jié)約計算成本。在文獻[21]中根據(jù)設(shè)定的典型系數(shù)（Canonical Coefficients）閾值選取SNP位點。當(dāng)有K個品種則需要K-1個閾值，隨著品種數(shù)的增多閾值的個數(shù)也會相應(yīng)增多。而本文都是由統(tǒng)一的排名方式選取SNP位點，而不需要特意設(shè)定閾值，對于數(shù)據(jù)的通用性更好。從表中可以看出在相同準確性的情況下，本文使用最少的SNP位點個數(shù)，以此提高計算速度，節(jié)約計算成本。

表2 不同方法的分類性能對比

4 結(jié)語

本文提出了利用PCA和隨機森林相結(jié)合的方法，可以從龐大的數(shù)據(jù)集中篩選出含有高信息量的SNP位點用于品種分類。實驗證實，對于將未知個體分配到與其對應(yīng)的品種的正確率達到了97.35%～98.37%。與其他SNP篩選方法相比較，將PCA和隨機森林很好地結(jié)合到一起，有效地減少用于分類的SNP位點數(shù)目，節(jié)約計算成本。并且本文方法即使隨機森林使用較少的SNP位點，在分類上也有良好的表現(xiàn)。但是對于其他物種和數(shù)據(jù)是否也有很好的性能還需要進一步的研究。