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

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ù)重點(diǎn)開放實(shí)驗(yàn)室，新疆維吾爾自治區(qū)動(dòng)物生物技術(shù)重點(diǎn)實(shí)驗(yàn)室，烏魯木齊 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)用價(jià)值。品種鑒別促進(jìn)了對遺傳信息的有效管理[1]；為育種策略的制定與實(shí)施奠定良好基礎(chǔ)[2]；為生物品牌產(chǎn)品的認(rèn)證提供了有效信息[3]；更為解決食品安全問題開辟了新途徑[4]。在以往品種鑒別研究中少有涉及羊的種類鑒別，而我國有著飼養(yǎng)羊群的悠久歷史，擁有豐富的羊群遺傳資源，據(jù)《中國畜禽遺傳資源志·羊志》記載的我國地方羊品種就達(dá)42個(gè)[5]，約占世界羊品種的5%[1]。本文將SNP位點(diǎn)應(yīng)用于羊的品種分類。

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

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

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

2 PCA和隨機(jī)森林相結(jié)合提取SNP位點(diǎn)

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

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

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

2.2 基于PCA的SNP位點(diǎn)提取

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

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

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

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

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

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

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

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

2.3 用于分類的隨機(jī)森林算法

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

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

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

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

2.4 位點(diǎn)重要性評價(jià)準(zhǔn)則

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

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

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

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

3 實(shí)驗(yàn)及結(jié)果分析

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

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

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

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

（3）訓(xùn)練集與測試集。從每個(gè)品種中隨機(jī)抽取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版本進(jìn)行PCA計(jì)算。為了減少計(jì)算時(shí)間，分別按照每一條染色體進(jìn)行計(jì)算，提高計(jì)算效率。每一個(gè)主成分所提取的信息量用方差來度量，其中某一個(gè)主成分方差的貢獻(xiàn)就等于原指標(biāo)相關(guān)矩陣相應(yīng)的特征值δi，則第i個(gè)主成分的方差貢獻(xiàn)率見式（6）：

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

圖3 所有SNP位點(diǎn)在前10個(gè)主成分上的方差貢獻(xiàn)率

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

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

圖4 520個(gè)SNP位點(diǎn)的前兩個(gè)主成分分析圖

圖5 520個(gè)SNP位點(diǎn)的前3個(gè)主成分分析圖

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

3.3 隨機(jī)森林重要性評價(jià)及分類模型建立

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

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

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

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

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

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

圖7 48個(gè)SNP位點(diǎn)的Gini指數(shù)值

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

圖8 被選中的SNP位點(diǎn)在每個(gè)染色體上的分布情況

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

為驗(yàn)證PCA可預(yù)選SNP位點(diǎn)，隨機(jī)森林可以進(jìn)一步減少用于分類的SNP位點(diǎn)個(gè)數(shù)并達(dá)到較好的準(zhǔn)確率?；?種數(shù)據(jù)集分別隨機(jī)抽取48和96個(gè)SNP位點(diǎn)對隨機(jī)森林的分類結(jié)果進(jìn)行測試，結(jié)果見表1。

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

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

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

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

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

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

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

4 結(jié)語

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