張小丹，李 喆，衛(wèi)澤剛*，劉 策，余凱哲，魏月華

（寶雞文理學(xué)院物理與光電技術(shù)學(xué)院，陜西 寶雞）

引言

隨著高通量測序技術(shù)的快速發(fā)展，產(chǎn)生的生物序列數(shù)據(jù)量呈指數(shù)級增加。面對海量序列數(shù)據(jù)，如何高效分析并挖掘隱藏在序列數(shù)據(jù)中的生物信息面臨新的挑戰(zhàn)[1]。其中序列比對是序列分析的前提和基礎(chǔ)，序列比對是按照一定規(guī)則，通過插入、刪除等操作進(jìn)行的有規(guī)律排列。通過序列比對，可以方便確定兩個或多個序列之間的相似性，進(jìn)而推斷序列間的同源性[2-3]。

Needleman-Wunsch（NW）算法是雙序列比對方法中最經(jīng)典的算法之一，由Saul Needlman 和Christian Wunsch 兩位科學(xué)家于1970 年提出[4]。借助于動態(tài)規(guī)劃打分策略，NW 算法可以記錄序列間的最優(yōu)比對得分，通過路徑回溯得到詳細(xì)的序列比對結(jié)果。在NW序列比對過程中，需要構(gòu)建M×N 大小的打分矩陣（M和N 分別是兩條序列的長度），其計(jì)算復(fù)雜度與空間復(fù)雜度均為O（MN），隨著序列長度的增加，具有很高的時間復(fù)雜度與空間復(fù)雜度[5-6]。因此，NW 算法不適合處理海量序列數(shù)據(jù)。

為快速計(jì)算序列間的相似性，研究者開發(fā)了許多基于序列k-mer 詞頻的“非比對”（alignment-free）算法[7]。k-mer 表示長度為k 的子片段，通過設(shè)定相應(yīng)的k 值，可以得到每條序列包含相應(yīng)子片段出現(xiàn)的次數(shù)（或稱為頻數(shù)），然后將序列轉(zhuǎn)化為k-mer 頻數(shù)（詞頻）向量，進(jìn)而借助其他向量相似性計(jì)算方法，得到序列間的k-mer 詞頻相似性。目前有很多種向量相似性計(jì)算方法，如歐氏距離、標(biāo)準(zhǔn)化歐氏距離、夾角余弦距離等。不同方法會得到不同的相似性結(jié)果，如何選取更接近于NW 算法的向量相似性方法，需要詳細(xì)的比較分析。因此，本文基于DNA 序列k-mer 詞頻向量，對目前常用的九種距離計(jì)算公式進(jìn)行總結(jié)，然后選取兩種數(shù)據(jù)集進(jìn)行測試，并與標(biāo)準(zhǔn)的NW 算法相似性進(jìn)行整體比較，分析不同距離計(jì)算方法的與標(biāo)準(zhǔn)NW 算法的差異，進(jìn)而幫助研究者選取合適的方法。

1 DNA 序列k-mer 詞頻向量表示

在運(yùn)用不同向量相似性計(jì)算方法之前，需要將DNA 序列轉(zhuǎn)成k-mer 詞頻向量。DNA 序列由A、C、G、T 四種堿基組成，因此，對于給定的k 值（k-mer 長度），k-mer 的可能性組合共有L 種（L=4k），然后計(jì)算每一種k-mer 在序列中出現(xiàn)的次數(shù)，即可構(gòu)建每條序列的k-mer 詞頻向量，向量長度即為L。假如現(xiàn)有兩條序列s 和t，其k-mer 詞頻向量可表示為hs=（sw1，sw2，…，swL）和ht=（tw1，tw2，…，twL），然后即可根據(jù)下面介紹的不同距離（相似性）計(jì)算方法得到hs與ht之間的k-mer 詞頻相似性。

2 常用距離計(jì)算方法

接下來介紹九種常用距離計(jì)算方法，這些方法大多數(shù)都可以在MATLAB 軟件[8-9]中找到相應(yīng)的計(jì)算函數(shù)并直接運(yùn)行，方便調(diào)用與計(jì)算。

2.1 歐氏距離

歐氏距離，又稱歐幾里得距離（Euclidean distance），是最常用、最易于理解的一種距離計(jì)算方法，源自歐式空間中兩點(diǎn)間的距離公式，衡量的是多維空間中兩個點(diǎn)之間的絕對距離，只要給出兩個數(shù)據(jù)的向量表示形式，即可計(jì)算得到歐氏距離。對于兩條序列s 和t，其k-mer 詞頻向量為hs=(sw1,sw2,...,swk)和ht=(sw1,sw2,...,swk)，序列間歐式距離的計(jì)算公式為：

2.2 標(biāo)準(zhǔn)化歐式距離

考慮到數(shù)據(jù)不同維度之間的尺度標(biāo)準(zhǔn)不一樣，標(biāo)準(zhǔn)化歐氏距離（Standard Euclidean distance）對上述歐式距離進(jìn)行了改進(jìn)，首先計(jì)算每個k-mer 分量的均值mean 和方差std，然后再用歐式距離進(jìn)行計(jì)算，即：

其中mean 是所有序列每個維度的平均值向量，std 是每個維度的標(biāo)準(zhǔn)差。

2.3 曼哈頓距離

曼哈頓距離（Manhattan distance）又稱為“折線距離”或“直角距離”，由十九世紀(jì)的德國數(shù)學(xué)家赫爾曼·閔可夫斯基提出，定義為兩個點(diǎn)在標(biāo)準(zhǔn)坐標(biāo)系上絕對軸距差總和，即每個分量間差值的絕對值之和。相對于歐式距離，計(jì)算更加簡單。對于兩條序列s 和t，其曼哈頓距離計(jì)算公式為：

2.4 切比雪夫距離

切比雪夫距離（Chebychev distance）可以認(rèn)為是曼哈頓距離的簡化版本，定義為兩個向量各個維度之間的最大距離差值。對于兩條序列s 和t，其切比雪夫距離計(jì)算公式為：

2.5 漢明距離

漢明距離（Hamming distance）一開始主要表示兩個等長字符串在每個位置對應(yīng)不同字符的數(shù)目，度量了通過替換字符的方式將其中一個字符串變成另一個字符串所需要的最小替換次數(shù)。對于本文中的序列詞頻向量，漢明距離通過判斷每個維度的數(shù)值是否相等進(jìn)行計(jì)算，定義為hs和ht中每個維度數(shù)值不相等的個數(shù)。

2.6 杰卡德距離

杰卡德距離（Jaccard distance）一開始用來衡量兩個集合之間差異性，定義為兩個集合的交集元素個數(shù)占并集元素個數(shù)的比例。對于詞頻向量hs和ht，杰卡德距離計(jì)算公式為：

2.7 余弦距離

幾何中夾角余弦可用來衡量兩個向量方向之間的差異，夾角余弦取值范圍為[-1，1]。夾角余弦越大表示兩個向量的夾角越小，夾角余弦越小表示兩向量的夾角越大。與以上的距離度量相比，夾角余弦更加注重兩個向量在方向上的差異。在機(jī)器學(xué)習(xí)或相似性度量中，從夾角余弦角度反映向量之間的差異性，可將兩個向量之間的夾角余弦值作為一種距離度量，稱為余弦距離（Cosine distance），計(jì)算公式為：

2.8 皮爾森相關(guān)系數(shù)

皮爾遜相關(guān)系數(shù)（Pearson correlation coefficient）主要用于度量兩個變量之間的相關(guān)程度，其值介于-1與1 之間。相關(guān)系數(shù)的絕對值越大，相關(guān)性越強(qiáng)。在機(jī)器學(xué)習(xí)或相似性度量中，皮爾遜相關(guān)系數(shù)也可作為一種距離方式，對于詞頻向量hs和ht，其皮爾遜相關(guān)距離計(jì)算公式為：

2.9 斯皮爾曼相關(guān)系數(shù)

斯皮爾曼相關(guān)系數(shù)（Spearman correlation coefficient）重點(diǎn)關(guān)注的是兩個變量之間單調(diào)關(guān)系的強(qiáng)度，即變化趨勢的強(qiáng)度。對于兩條序列的k-mer 向量hs和ht，首先需要對hs和ht進(jìn)行排序，然后得到每個向量的秩（次序）向量ds和dt，最后根據(jù)下面公式計(jì)算得到斯皮爾曼相關(guān)距離：

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

3.1 評價指標(biāo)

本文首先計(jì)算每個方法的距離，然后轉(zhuǎn)換成相似性。對于相似性大于1 的方法，需要對其進(jìn)行標(biāo)準(zhǔn)化，即每個計(jì)算結(jié)果減去最小值，再除以最大值與最小值之差。這樣可以將每個方法的相似性歸一化到0 到1之間，方便比較。然后通過繪制不同序列相似性計(jì)算方法與標(biāo)準(zhǔn)NW 相似性的散點(diǎn)圖來直觀比較每個方法與標(biāo)準(zhǔn)NW 方法相似性的相關(guān)程度。然后添加擬合直線，可以進(jìn)一步量化相關(guān)性強(qiáng)度，幫助分析不同非比對方法與標(biāo)準(zhǔn)NW 相似性之間的符合程度。最后測試了每個方法的運(yùn)行時間，比較不同方法計(jì)算時間的快慢。

3.2 K 值選取

K 值（k-mer 長度）選取會影響序列k-mer 詞頻向量的長度，進(jìn)而影響最終的相似性計(jì)算。如果設(shè)置較小的k 值，相應(yīng)的詞頻向量長度也較短，不能提取足夠的序列信息。如果設(shè)置較大的k 值，會增加向量長度，降低計(jì)算效率。因此，需要根據(jù)數(shù)據(jù)集選取合適的k 值，文獻(xiàn)[10]給出了k 值的計(jì)算公式[10]，如下所示：

式中，n 為數(shù)據(jù)集S 中的序列的個數(shù)，len(i)表示第i 條序列的長度，ceil 表示向上取整。因此，對于給定的數(shù)據(jù)集，先通過公式計(jì)算k 值，然后再將每條序列轉(zhuǎn)成相應(yīng)的k-mer 詞頻向量。

為了比較不同距離計(jì)算方法與NW 算法間的差異，本文選取兩組數(shù)據(jù)進(jìn)行比較分析：模擬數(shù)據(jù)集和病毒數(shù)據(jù)集。

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

首先采用一組DNA 模擬數(shù)據(jù)集進(jìn)行比較分析[11]，是序列聚類與比對中常用的測試數(shù)據(jù)，總共包含236條序列，序列長度范圍為998～1 037，平均序列長度為1 016。

首先根據(jù)公式計(jì)算適用于此數(shù)據(jù)的k 值，計(jì)算結(jié)果為k=4，然后將每一條序列轉(zhuǎn)換為相應(yīng)的k-mer 詞頻向量，最后根據(jù)不同的計(jì)算公式得到相似性計(jì)算結(jié)果。每個方法的相似性計(jì)算結(jié)果與標(biāo)準(zhǔn)NW 相似性的散點(diǎn)圖如圖1 所示，其中橫坐標(biāo)表示標(biāo)準(zhǔn)NW 相似性，縱坐標(biāo)是不同方法的相似性，每個圖中的直線是相應(yīng)一次擬合直線，可通過MATLAB 相關(guān)函數(shù)計(jì)算得到。從圖1 中可以觀察到，整體來看，相對于其他方法，歐氏距離、標(biāo)準(zhǔn)歐氏距離、曼哈頓距離和余弦距離的直線擬合相關(guān)系數(shù)更接近于1，說明這四個方法的距離計(jì)算結(jié)果與標(biāo)準(zhǔn)NW 算法更接近。其他方法的計(jì)算結(jié)果與NW 相似性的相關(guān)系數(shù)相差0.3 左右，與NW 相似性整體偏差較大。

圖1 基于k-mer 詞頻的不同方法針對模擬數(shù)據(jù)序列的相似性計(jì)算結(jié)果

接下來用病毒數(shù)據(jù)測試每個方法對長序列的計(jì)算結(jié)果，此數(shù)據(jù)集同樣來自文獻(xiàn)[10]，共包含96 條序列，序列長度范圍為2 605～13 246，平均序列長度為6 624，可以測試不同方法對于長序列數(shù)據(jù)的計(jì)算結(jié)果。每個方法與標(biāo)準(zhǔn)NW 距離的散點(diǎn)圖分布如圖2 所示?？梢钥闯觯瑲W氏距離與杰拉德距離的整體相關(guān)程度最接近于1，說明這兩個方法在處理長病毒序列時與NW 更接近，其中杰拉德的分布更集中在擬合直線附近，而歐式距離分布較分散，說明杰拉德距離在計(jì)算長病毒序列比歐式距離更好。

圖2 基于k-mer 詞頻的不同方法針對病毒數(shù)據(jù)集序列相似性計(jì)算結(jié)果

表1 是不同方法針對模擬數(shù)據(jù)和病毒數(shù)據(jù)的計(jì)算時間，可以看出，針對模擬數(shù)據(jù)，NW 算法計(jì)算每對序列間的相似性需要723 秒，而其他方法用時均不到0.1 秒，針對病毒數(shù)據(jù)，NW 算法需要6 178 秒，而其他方法用時均不到2 秒。說明基于k-mer 詞頻的計(jì)算方法至少比NW 算法快103倍，因此，基于k-mer 詞頻方法的計(jì)算時間更加高效，可以極大降低計(jì)算時間。

表1 不同方法計(jì)算時間

4 結(jié)論

生物序列比對是生物信息學(xué)研究領(lǐng)域的基礎(chǔ)分析工作，可為序列相似性分析奠定重要的理論依據(jù)。傳統(tǒng)的序列相似性分析方法通常是基于比對的算法，具有很高的時間和空間復(fù)雜度高。隨著測序技術(shù)的快速發(fā)展，生物序列以指數(shù)級的形式快速增長，面對目前海量級的生物數(shù)據(jù)，標(biāo)準(zhǔn)雙序列比對NW 算法在計(jì)算序列相似性時具有較高的時間復(fù)雜度。為快速得到序列間的相似性，需要借助于k-mer 詞頻向量方法。本文從序列k-mer 詞頻向量及常用的九中k-mer 詞頻計(jì)算方法進(jìn)行了詳細(xì)介紹，并用兩種數(shù)據(jù)集進(jìn)行了比較分析。實(shí)驗(yàn)結(jié)果表明，基于k-mer 詞頻相似性計(jì)算方法比標(biāo)準(zhǔn)NW 算法速度至少快103倍，但不同的k-mer 詞頻計(jì)算方法得到的相似性與標(biāo)準(zhǔn)NW 算法差別較大，相對而言，歐式距離在兩個數(shù)據(jù)集的相似性結(jié)果與NW 更接近，在計(jì)算大規(guī)模序列相似性時，可以作為優(yōu)先選擇的方法。