基于k- mer 詞頻向量的九種DNA 序列相似性計算方法比較分析

張小丹，李 喆，衛澤剛*，劉 策，余凱哲，魏月華

（寶雞文理學院物理與光電技術學院，陜西 寶雞）

引言

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

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

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

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

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

2 常用距離計算方法

接下來介紹九種常用距離計算方法，這些方法大多數都可以在MATLAB 軟件[8-9]中找到相應的計算函數并直接運行，方便調用與計算。

2.1 歐氏距離

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

2.2 標準化歐式距離

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

其中mean 是所有序列每個維度的平均值向量，std 是每個維度的標準差。

2.3 曼哈頓距離

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

2.4 切比雪夫距離

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

2.5 漢明距離

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

2.6 杰卡德距離

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

2.7 余弦距離

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

2.8 皮爾森相關系數

皮爾遜相關系數（Pearson correlation coefficient）主要用于度量兩個變量之間的相關程度，其值介于-1與1 之間。相關系數的絕對值越大，相關性越強。在機器學習或相似性度量中，皮爾遜相關系數也可作為一種距離方式，對于詞頻向量hs和ht，其皮爾遜相關距離計算公式為：

2.9 斯皮爾曼相關系數

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

3 實驗結果

3.1 評價指標

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

3.2 K 值選取

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

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

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

3.3 實驗結果

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

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

圖1 基于k-mer 詞頻的不同方法針對模擬數據序列的相似性計算結果

接下來用病毒數據測試每個方法對長序列的計算結果，此數據集同樣來自文獻[10]，共包含96 條序列，序列長度范圍為2 605～13 246，平均序列長度為6 624，可以測試不同方法對于長序列數據的計算結果。每個方法與標準NW 距離的散點圖分布如圖2 所示。可以看出，歐氏距離與杰拉德距離的整體相關程度最接近于1，說明這兩個方法在處理長病毒序列時與NW 更接近，其中杰拉德的分布更集中在擬合直線附近，而歐式距離分布較分散，說明杰拉德距離在計算長病毒序列比歐式距離更好。

圖2 基于k-mer 詞頻的不同方法針對病毒數據集序列相似性計算結果

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

表1 不同方法計算時間

4 結論

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