基于位置信息的DNA 序列特征提取＊

深圳大學(xué) 陳煜元 周小安

DNA 序列的分類是生物信息學(xué)的主要研究任務(wù)之一，如何提取DNA 序列中的特征是影響分類精度的重要因素。為了更好地保留序列中堿基的信息，本文提出了一種基于堿基距離和相關(guān)性的特征提取方法。以H1N1、H5N1、COVID-19 等6 種病毒作為研究對(duì)象，將DNA序列轉(zhuǎn)化為特征向量，并用KNN 算法對(duì)冠狀和非冠狀病毒進(jìn)行分類。實(shí)驗(yàn)結(jié)果表明該方法能提高分類的準(zhǔn)確率。

據(jù)估計(jì)地球上約有1000 萬(wàn)～1 億種生物，如此龐大的數(shù)據(jù)使得生物分類面臨著巨大挑戰(zhàn)[1]，因此DNA 序列的分類成為了人們的研究熱點(diǎn)，也是當(dāng)前生物信息學(xué)的主要研究任務(wù)之一。特征提取是DNA 序列分類研究中至關(guān)重要的一環(huán)，旨在最大限度保留原序列數(shù)據(jù)的基礎(chǔ)上將序列轉(zhuǎn)化為數(shù)值特征，以挖掘其中所存在的生物規(guī)律。隨著計(jì)算機(jī)技術(shù)的發(fā)展和測(cè)序技術(shù)的不斷進(jìn)步，堿基的組成和分布信息在DNA 序列特征提取中備受關(guān)注[2]。最基本的特征提取方法為K-mers[3]，該方法隨著k 的增大特征維數(shù)呈現(xiàn)指數(shù)級(jí)的增長(zhǎng)，而在訓(xùn)練樣本不足的情況下高維數(shù)據(jù)的研究會(huì)帶來(lái)“過(guò)擬合”“維數(shù)災(zāi)難”等問題[4]，故k 的取值不能太大，而特征維數(shù)不足可能會(huì)丟失序列中的重要信息。此外，K-mers 方法忽略了堿基的距離和排列情況[5]。因此，本文擬提取出基于相同堿基間距離和不同堿基間相關(guān)性的特征用于病毒序列分類。該特征提取方法以DNA 序列中堿基的位置為基礎(chǔ)，分別記錄各堿基出現(xiàn)的位置，再通過(guò)合適的數(shù)學(xué)方法計(jì)算出平均距離和相關(guān)系數(shù)。實(shí)驗(yàn)結(jié)果表明，新的特征提取方法在KNN 分類器上能取得較好的分類效果。

1 KNN 算法

1.1 KNN 簡(jiǎn)介

K 近鄰（K-Nearest Neighbor）算法簡(jiǎn)稱KNN，是Cover 和Hart 在1968 年時(shí)首先提出的。它是一個(gè)在理論上較為成熟的算法，也是最常用、最簡(jiǎn)單的機(jī)器學(xué)習(xí)算法之一。由于和其他分類算法相比沒有顯示的學(xué)習(xí)過(guò)程，所依據(jù)的“多數(shù)決定”的思想很容易理解，在多分類問題上表現(xiàn)的比其他分類算法要好，而且計(jì)算過(guò)程經(jīng)過(guò)優(yōu)化后能夠大幅降低計(jì)算次數(shù)，因此在分類領(lǐng)域有著廣泛的應(yīng)用[6]。算法的原理是將待分類的樣本與訓(xùn)練集的樣本逐一計(jì)算出距離，按距離從小到大進(jìn)行排序，然后取出最近的K 個(gè)訓(xùn)練樣本，這K 個(gè)樣本中數(shù)量較多那一類即為測(cè)試樣本的類別。

1.2 距離的計(jì)算

KNN 計(jì)算序列樣本之間距離的方法主要有曼哈頓距離、歐式距離和閔可夫斯基距離等，本文采用的是歐式距離法。令訓(xùn)練樣本為X，特征的向量表示為(x1,x2,…,xn)，測(cè)試樣本為Y，特征為(y1,y2,…,yn)，則它們的歐式距離如式（1）所示：

式中n 表示X 和Y 的特征維數(shù)，若n=2 相當(dāng)于計(jì)算平面上兩點(diǎn)之間的距離，n=3 相當(dāng)于計(jì)算三維空間中兩點(diǎn)的距離。上式計(jì)算的是一個(gè)訓(xùn)練樣本和一個(gè)測(cè)試樣本之間的距離，實(shí)際的分類問題中往往有m 個(gè)訓(xùn)練樣本{X1,X2,…,Xm}，需分別計(jì)算待分類樣本與每個(gè)訓(xùn)練樣本的距離d(Xi,Y)(i=1,2,…,m)，從小到大排序后再進(jìn)行下一步工作。

1.3 K 值的選取

K 值的選取是KNN 算法中至關(guān)重要的一環(huán)，取值太小容易導(dǎo)致過(guò)擬合，太大則會(huì)使得分類誤差增大。不僅如此，K 值有時(shí)能直接影響到分類結(jié)果。如圖1 所示，假設(shè)一個(gè)二分類問題，藍(lán)色正方形代表類別A，綠色三角形代表類別B，紅色圓形為待分類樣本X。若取K=3，即把距離最近的3 個(gè)樣本作為依據(jù)，此時(shí)類別B 的個(gè)數(shù)為2，類別A 的個(gè)數(shù)為1，此時(shí)KNN 會(huì)將待分類樣本歸為類別B；而取K=5 時(shí)顯然類別A 的個(gè)數(shù)比類別B 的多，則X 會(huì)被歸為A 類。

圖1 K 近鄰算法Fig.1 K-nearest neighbors

為了降低上述情況帶來(lái)的影響，可以采用距離加權(quán)的方式，給每個(gè)已知類別的樣本賦予權(quán)重，其值和訓(xùn)練樣本與測(cè)試樣本之間距離成反比，這樣就使得較相似的樣本點(diǎn)在分類上具有更高的權(quán)重。本文將通過(guò)經(jīng)驗(yàn)法選取合適的K值，即對(duì)不同的K 值進(jìn)行重復(fù)實(shí)驗(yàn)，得到一個(gè)使分類準(zhǔn)確率最高的結(jié)果。

2 序列特征的提取

2.1 基于堿基頻率的特征提取

這是一種比較典型的特征提取方法，即k-mers法，以各堿基或堿基組合在DNA 序列中的頻率作為特征。由于堿基的種類有4種，故對(duì)于k(k=1,2,3,…)個(gè)堿基的組合，共有4k種組合方式。

單堿基有A、T、G、C 4種，在不同DNA 序列甚至同一序列的不同片段中，每種堿基出現(xiàn)的頻率也是不相同的。設(shè)一條含有m 個(gè)堿基的序列，其中堿基n 的個(gè)數(shù)為c，則堿基n 的頻率可記為Pn=c/m。把4 種堿基的頻率用向量表示如式（2）所示：

其中Pi(i∈{A,T,G,C})表示堿基i 在序列中出現(xiàn)的頻率。

雙堿基有AA、AT、AG、AC、TA、TT、TG、TC、GA、GT、GG、GC、CA、CT、CG、CC 共42種，即16種組合方式。實(shí)驗(yàn)采用滑動(dòng)窗口算法對(duì)各雙堿基的個(gè)數(shù)進(jìn)行計(jì)算，這樣可以降低堿基缺失對(duì)實(shí)驗(yàn)結(jié)果的影響。一條含有m 個(gè)堿基的序列，用滑動(dòng)窗口法得到的雙堿基共有m-1個(gè)，故雙堿基n 的頻率可記為Pn=c/(m-1)，c為雙堿基n 出現(xiàn)的次數(shù)。可將雙堿基在序列中的頻率表示為如下16 維向量：

其中Pij(i,j∈{A,T,G,C})表示堿基ij 在序列中出現(xiàn)的頻率。

基于三堿基的表示方法共有64 種組合，同樣用滑動(dòng)窗口法計(jì)算出各三堿基的頻率，將其用如下64 維向量表示：

其中Pijk(i,j,k∈{A,T,G,C})表示堿基ijk 在序列中出現(xiàn)的頻率。

基于四堿基組合的表示方法共有44即256種，五堿基的共有1024種，隨著n 的增大組合方式呈指數(shù)型增長(zhǎng)。在實(shí)際的機(jī)器學(xué)習(xí)算法中，太多的特征會(huì)使得計(jì)算的時(shí)間成本增加，且可能導(dǎo)致過(guò)擬合。此外，多個(gè)堿基的組合在序列中出現(xiàn)頻率較低，故不考慮3 個(gè)堿基以上的頻率作為特征輸入。

上述基于堿基頻率的特征提取方法雖然能取得較好的分類效果，但是并不能表達(dá)出堿基的位置信息。本文提出一種新的DNA 序列特征提取方法，與k-mers 方法不同的是，新方法在考慮了堿基頻率的基礎(chǔ)上還包含了序列中堿基的距離和相關(guān)性信息。

2.2 基于堿基位置信息的特征提取

2.2.1 基于堿基距離的特征提取

僅用堿基的頻率還不足以描述一條DNA 序列，因?yàn)閮蓷l完全不同的序列可能會(huì)出現(xiàn)相同的堿基頻率，如以下兩條長(zhǎng)度為20 個(gè)堿基的序列：

序列a：ATCGC GCAGA GATAT CTATA

序列b：GCACA TCAGA TCAGA TATGT

序列a 和序列b 的A、T、G、C 含量完全相同，雙堿基含量和三堿基含量如AT、AG、TC、GCA、TAT等也有相同或相似之處，但這卻是截然不同的兩條序列。因此，為了使得分類更加準(zhǔn)確，可用一種新的基于堿基距離的特征表示方法，即堿基（或堿基組合）之間的平均距離來(lái)描述一條序列。設(shè)有一條長(zhǎng)度為m 的DNA 序列S=s1,s2,...sm，用Cn表示堿基n(n∈{A,T,G,C})的個(gè)數(shù)，Lni來(lái)表示堿基n 在序列中第i 次(i=1,2,…,Cn)出現(xiàn)的位置，則相鄰堿基n 之間的距離Dnj=Lnj+1-Lnj(j=1,2,…,Cn-1)，序列中堿基n 的平均距離為：

以序列a為例，堿基A的個(gè)數(shù)CA=7，位置分別為L(zhǎng)A1=1，LA2=8，LA3=10，LA4=12，LA5=14，LA6=18，LA7=20，相鄰兩個(gè)堿基A 之間的距離DA1=7，DA2=2，DA3=2，DA4=2，DA5=4，DA6=2，可算出其平均距離：

同理堿基T 分別位于第2，13，15，17，19處，故堿基T 的平均距離DT=3.4；用該方法算出堿基G 和堿基C的平均距離分別為DG=2.33，DC=4.33。對(duì)于序列b 也用此法計(jì)算得：DA’=2.3，DT’=3.5，DG’=6，DC’=3.33，即：

其中Da、Db分別表示序列a 和序列b 中堿基A、T、G、C 的平均距離向量。可看出在單堿基含量相同的情況下，堿基之間的平均距離可能會(huì)有較大的差異。

對(duì)雙堿基的平均距離和單堿基的計(jì)算方法類似，堿基組合的位置以第一個(gè)堿基為準(zhǔn)。例如要計(jì)算序列a 中堿基組合AT 的平均距離，可找出4 個(gè)AT，其中堿基A 的位置分別為1、12、14、18，故AT 的平均距離DAT=5.67；同理序列b 中AT 的平均距離DAT’=4。堿基AT 在序列a 和序列b 中含量相同（都是4 個(gè)），平均距離卻有所差異，這正說(shuō)明了我們不能只關(guān)注序列中堿基的含量而忽略了堿基的距離信息。

2.2.2 基于堿基相關(guān)性的特征提取

序列中不同堿基之間的相關(guān)關(guān)系也是區(qū)分不同種類生物的重要特征之一，將堿基n 的位置分布記作Sn(x)(x=1,2,…,m)，其定義為：

式中Lni表示堿基n 在DNA 序列中第i 次出現(xiàn)的位置，Cn表示堿基n 的個(gè)數(shù)。在堿基n 第一次出現(xiàn)之前，位置記為0；兩個(gè)堿基n 之間的位置記為它們的平均值；最后一個(gè)堿基n 之后的位置都記為最后一個(gè)堿基n 出現(xiàn)的位置。

以序列a 為例，CT=5，LT1=2，LT2=13，LT3=15，LT4=17，LT5=19，故ST(2)=2，ST(13)=13，ST(15)=15，ST(17)=17，ST(19)=19；第一個(gè)T 出現(xiàn)之前，ST(1)=0；最后一個(gè)T 之后，ST(20)=19；兩個(gè)T 之間，ST(3)=ST(4)=…=ST(12)=(2+13)/2=7.5，ST(14)=(13+15)/2=14，ST(16)=(15+17)/2=16，ST(18)=(17+19)/2=18，即堿基T 在序列a 中的位置分布為：

同理計(jì)算出堿基G 在序列a 中的位置分布為：

為體現(xiàn)序列中兩種不同堿基的相關(guān)程度并用數(shù)值表示出來(lái)，可用皮爾森相關(guān)系數(shù)來(lái)計(jì)算，定義是兩個(gè)變量之間的協(xié)方差和標(biāo)準(zhǔn)差的商：

綜上，本文提取出DNA 序列中基于堿基頻率、堿基間的距離和相關(guān)性的110 維特征向量：

將其作為KNN 算法中的特征輸入，對(duì)病毒序列進(jìn)行分類研究。

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

本文采用的實(shí)驗(yàn)數(shù)據(jù)均來(lái)源于美國(guó)國(guó)家生物技術(shù)信息中心（NCBI）。從NCBI 網(wǎng)站上下載6 種不同病毒(H1N1，H5N1，H7N9，SARS，MERS，COVID-19)的DNA 序列，前3 種為非冠狀病毒，后3 種為冠狀病毒。每種病毒各取50 組序列，總共300組，其中150 組的類別為“冠狀”，另外150 組類別為“非冠狀”，每組序列包含240 個(gè)堿基。用Python 中的scikit-learn 模塊對(duì)序列進(jìn)行二分類，評(píng)估指標(biāo)為分類準(zhǔn)確率（Accuracy），即正確分類的序列個(gè)數(shù)占總序列個(gè)數(shù)的比例。

首先用KNN 算法進(jìn)行兩輪實(shí)驗(yàn)，第一輪實(shí)驗(yàn)的數(shù)據(jù)是堿基的頻率特征，即序列中單堿基、雙堿基和三堿基的頻率，共84 維數(shù)據(jù)。第二輪實(shí)驗(yàn)的數(shù)據(jù)是在第一輪實(shí)驗(yàn)的基礎(chǔ)上增加堿基的距離和相關(guān)性特征，共110 維數(shù)據(jù)。兩輪實(shí)驗(yàn)分別測(cè)試不同K 值下的分類準(zhǔn)確率，結(jié)果如表1 所示。

表1 兩種特征提取方法的分類準(zhǔn)確率Tab.1 Classification accuracy of the two feature extraction methods

由表1 可知，在K 值為3～7 的情況下，分類準(zhǔn)確率呈拋物線趨勢(shì)，先增后減。在僅采用堿基頻率特征的算法中，K=5 時(shí)分類準(zhǔn)確率最高，為96.02%，K=4 次之；增加了堿基距離和相關(guān)性特征的實(shí)驗(yàn)中，K=6 時(shí)分類準(zhǔn)確率最高，為97.72%，K=5 次之。無(wú)論K 的取值為多大，增加了距離特征之后的分類模型都是更有效的，準(zhǔn)確率均提高1%～2%左右。

上述實(shí)驗(yàn)中訓(xùn)練樣本與測(cè)試樣本的比例為7:3，即210 組訓(xùn)練序列，90 組測(cè)試序列。考慮到測(cè)試樣本的歸類是以訓(xùn)練樣本的類別為依據(jù)的，此比例可能會(huì)影響分類結(jié)果。本文在訓(xùn)練樣本數(shù)量應(yīng)大于測(cè)試樣本數(shù)量的原則上，適當(dāng)調(diào)整了比例，以堿基的頻率和位置信息的110維特征作為輸入進(jìn)行了重復(fù)實(shí)驗(yàn)，結(jié)果如表2 所示。

由表2 可知，訓(xùn)練集與測(cè)試集比例為8:2 且K=6時(shí)，訓(xùn)練集與測(cè)試集比例為9:1 且K=5 或K=6時(shí)，分類的準(zhǔn)確率能達(dá)到98%以上。比例為8:2 且K 的取值為5 時(shí)的準(zhǔn)確率非常接近98%，這正說(shuō)明K 的取值偏大或偏小都會(huì)使得分類效果降低，在K=5 和K=6 時(shí)可達(dá)到較高的準(zhǔn)確率。此外，無(wú)論K 的取值為多少，分類準(zhǔn)確率都會(huì)隨著訓(xùn)練集和測(cè)試集比例的增加而提高，即訓(xùn)練樣本較多時(shí)模型能更加高效地學(xué)習(xí)。綜上，將采用訓(xùn)練集與測(cè)試集比例為9:1，K=6 的KNN 模型進(jìn)行序列的分類，此時(shí)準(zhǔn)確率最高，為98.47%。

表2 不同比例下的分類準(zhǔn)確率Tab.2 Classification accuracy of different proportions

4 結(jié)語(yǔ)

本文闡述了DNA 序列中堿基位置信息的重要性，提出一種基于堿基之間的距離和相關(guān)性特征表示方法，將其運(yùn)用于6 種病毒序列的特征提取，并利用KNN 算法進(jìn)行分類，實(shí)驗(yàn)結(jié)果表明該特征提取方法能提高分類準(zhǔn)確率。這種方法提取的特征向量維數(shù)較高，因此更適用于多個(gè)DNA 序列的分類研究。

引用

[1] 竇向梅,肖暉,黃大衛(wèi).DNA分類概述[J].生物學(xué)通報(bào),2008,43(6):23-25.

[2] 劉福樂.DNA、RNA和蛋白質(zhì)序列特征提取方法研究及應(yīng)用[D].哈爾濱:哈爾濱工業(yè)大學(xué),2015.

[3] Shobhit Gupta,Jonathan Dennis,Robert E Thurman,et al.Predicting Human Nucleosome Occupancy from Primary Sequence[J].PLoS Computational Biology,2008,4(8): e1000134.

[4] 李郅琴.特征選擇方法綜述[J].計(jì)算機(jī)工程與應(yīng)用,2019,55 (24):10-11.

[5] 韓軼平,余杭,劉威,等.DNA序列的分類[J].數(shù)學(xué)的實(shí)踐與汄識(shí),2001(1):38-45.

[6] 皮亞宸.K近鄰分類算法的應(yīng)用研究[J].通訊世界,2019,26(1): 286-287.