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

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

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

引言

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

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

為快速計(jì)算序列間的相似性，研究者開發(fā)了許多基于序列k-mer 詞頻的“非比對(duì)”（alignment-free）算法[7]。k-mer 表示長(zhǎng)度為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)化歐氏距離、夾角余弦距離等。不同方法會(huì)得到不同的相似性結(jié)果，如何選取更接近于NW 算法的向量相似性方法，需要詳細(xì)的比較分析。因此，本文基于DNA 序列k-mer 詞頻向量，對(duì)目前常用的九種距離計(jì)算公式進(jìn)行總結(jié)，然后選取兩種數(shù)據(jù)集進(jìn)行測(cè)試，并與標(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 四種堿基組成，因此，對(duì)于給定的k 值（k-mer 長(zhǎng)度），k-mer 的可能性組合共有L 種（L=4k），然后計(jì)算每一種k-mer 在序列中出現(xiàn)的次數(shù)，即可構(gòu)建每條序列的k-mer 詞頻向量，向量長(zhǎng)度即為L(zhǎng)。假如現(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)間的距離公式，衡量的是多維空間中兩個(gè)點(diǎn)之間的絕對(duì)距離，只要給出兩個(gè)數(shù)據(jù)的向量表示形式，即可計(jì)算得到歐氏距離。對(duì)于兩條序列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）對(duì)上述歐式距離進(jìn)行了改進(jìn)，首先計(jì)算每個(gè)k-mer 分量的均值mean 和方差std，然后再用歐式距離進(jìn)行計(jì)算，即：

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

2.3 曼哈頓距離

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

2.4 切比雪夫距離

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

2.5 漢明距離

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

2.6 杰卡德距離

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

2.7 余弦距離

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

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

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

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

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

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

3.1 評(píng)價(jià)指標(biāo)

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

3.2 K 值選取

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

式中，n 為數(shù)據(jù)集S 中的序列的個(gè)數(shù)，len(i)表示第i 條序列的長(zhǎng)度，ceil 表示向上取整。因此，對(duì)于給定的數(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]，是序列聚類與比對(duì)中常用的測(cè)試數(shù)據(jù)，總共包含236條序列，序列長(zhǎng)度范圍為998～1 037，平均序列長(zhǎng)度為1 016。

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

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

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

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

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

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

4 結(jié)論

生物序列比對(duì)是生物信息學(xué)研究領(lǐng)域的基礎(chǔ)分析工作，可為序列相似性分析奠定重要的理論依據(jù)。傳統(tǒng)的序列相似性分析方法通常是基于比對(duì)的算法，具有很高的時(shí)間和空間復(fù)雜度高。隨著測(cè)序技術(shù)的快速發(fā)展，生物序列以指數(shù)級(jí)的形式快速增長(zhǎng)，面對(duì)目前海量級(jí)的生物數(shù)據(jù)，標(biāo)準(zhǔn)雙序列比對(duì)NW 算法在計(jì)算序列相似性時(shí)具有較高的時(shí)間復(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 算法差別較大，相對(duì)而言，歐式距離在兩個(gè)數(shù)據(jù)集的相似性結(jié)果與NW 更接近，在計(jì)算大規(guī)模序列相似性時(shí)，可以作為優(yōu)先選擇的方法。