基于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 更接近，在计算大规模序列相似性时，可以作为优先选择的方法。