基于前缀标识符及其位置的DNA序列比较

Comparison of DNA Sequences Based on Prefix Identifiers and Their Locations

摘要: 分子序列比较是生物信息学中最基本、最主要的问题，DNA序列相似性分析是研究的重要的课题。非比对方法是研究序列比较的方法之一，它克服了比对方法的局限，其计算速度更快。本文从前缀标识符位置角度出发，利用信息熵，提出了序列分析的非比对方法。本文通过对生物序列构建前缀树，得到生物序列前缀标识符的基础上，以两两序列的共同前缀标识符为研究对象，提取它们在序列中位置信息，将它们的位置差的绝对值看成随机变量，利用信息熵，提出新的DNA序列相似性度量方法，建立有效的模型。将70个哺乳动物的线粒体DNA序列作为实验数据集，应用该模型得到的相似性距离构建生物进化树。该进化树的分类结果符合当前的生物学分类标准。

Abstract: Comparison of molecular sequence is the most basic and important problem in bioinformatics. DNA sequence similarity analysis is an important research topic. Alignment-free method is one of the methods to study sequence comparison. It overcomes the limitation of alignment method and is faster than alignment method. In this paper, from the point of view of prefix identifier location, the alignment-free method of sequence analysis is proposed by using information entropy. Based on the prefix tree and the prefix identifier of biological sequences, the position information of pairwise sequences is extracted by using the common prefix identifiers of pairwise sequences. The absolute value of their position difference is regarded as random variable. Using information entropy, a new DNA sequence similarity measurement method is proposed and an effective model is established. Mitochondrial DNA sequences of 70 mammalian were used as experimental data sets. Construct the Phylogenetic tree based on the similarity distance obtained by the model. The classification results of Phylogenetic tree conform to the current biological classification.

1. 引言

近几十年来，随着人类基因组计划的实施和完成，分子生物学发展呈现出生物数据爆炸式增长的特点。面对呈指数增长趋势的海量数据，高效管理、准确解读、从而挖掘有用的生物信息是一项有意义的工作，同时也是生物、数学、计算机科学等多个领域专家学者面临的一大挑战 [1]。本文将信息理论知识、前缀树知识与生物信息学知识相结合，提出一种新的序列相似性度量方法。

1973年Peter Weiner提出了前缀树模型及其算法 [2]。前缀树模型的提出启发很多在生物信息学领域的学者们，为非比对方法提供了新的研究手段。部分文献基于前缀树，利用序列最长公共子串，提出序列相似性度量模型 [3]。有文献通过讨论k词在序列中的位置分布，来研究序列相似性 [4]。本文将提出一种新的以DNA序列前缀标识符为研究对象的非比对方法。

本文将提出的新的DNA序列相似性度量模型在70条哺乳动物线粒体DNA序列上进行实验，根据生成的进化树，进行结果分析与讨论。将所得的结果首先与NCBI给出的生物学分类标准进行比较，证明本文提出的方法对该数据集具有正确的分类能力，其次与已发表文献结果进行比较，得出本结果在不同程度上由于前人结果。

2. 材料

本文中使用哺乳动物线粒体DNA序列作为数据集进行实验。本数据集中共有70条不同的哺乳动物线粒体DNA序列。序列长度最短为16295，最长为17447。并且DNA序列中不含有非A、C、G、T的字符。本文中的数据集均为从网站National Center for Biotechnology Information“NCBI”(https://www.ncbi.nlm.nih.gov/)中Genbank数据库下载得到。具体信息见附录。

本数据集中序列来自于哺乳动物中原兽亚纲动物(Prototheria)、有袋类动物(Marsupialia)、有胎盘类动物(Placentalia)。其中有胎盘类哺乳动物的序列来自于非洲总兽目(Afrotheria)、贫齿目(Xenarthra)、灵长总目(Euarchontoglires)、劳亚兽总目(Laurasiatheria)四个超目。本数据集中有17个目，41个科。

本文数据集与文献 [5] 有63个相同物种，有3条序列有相同的属，2条序列有相同的科。

3. 方法

3.1. 前缀标识符

3.1.1. 线形序列前缀标识符定义

本文对前缀标识符的定义来源于文献 [1]。对任意序列S， $S=S\left(1\right)S\left(2\right)\cdots S\left(n\right)$，则n为序列S的长度，其中 $S\left(i\right)$ 称为序列S在位置i处的字符。将字符串 $S\left(i\right)S\left(i+1\right)S\left(i+2\right)\cdots S\left(j\right)$ 称为序列S中i处开始在j处结束的字符串。在序列S的结尾处添加一个结束字符，该字符必须与S中每一个字符都不相同，例如#。于是得到连接后的序列，记为SS，则 $SS=S\left(1\right)S\left(2\right)\cdots S\left(n\right)\#$。若字符串 $S\left(i\right)S\left(i+1\right)\cdots S\left(j\right)$ 为位置i处的前缀标识符，则满足以下两个条件：(1) 字符串 $S\left(i\right)S\left(i+1\right)\cdots S\left(j\right)$ 序列SS中只出现一次(2) $S\left(i\right)S\left(i+1\right)\cdots S\left(j-1\right)$ 在序列SS中至少出现两次。例如，字符串 $S=acppb$，添加结束字符#，则 $SS=acppb\#$。那么S在位置1处的前缀标识符为a，位置3处的前缀标识符为pp，而不是p，这是因为在SS中位置3和位置4处均出现字符p，即不满足条件(1)，pp在SS只出现一次，而p在SS中出现2次，满足条件(1)和(2)。

3.1.2. 环形序列前缀标识符定义

因为哺乳动物线粒体DNA序列是环形的，因此本文定义环形序列的前缀标识符，其定义如下：任意DNA序列S1， $S1=S\left(1\right)S\left(2\right)\cdots S\left(n\right)$，若字符串 $S\left(i\right)S\left(i+1\right)\cdots S\left(j\right)$ 为位置i处的前缀标识符，则满足以下两个条件：(1) 在 $S\left(i\right)S\left(i+1\right)\cdots S\left(j\right)$ 序列S1中只出现一次(2) $S\left(i\right)S\left(i+1\right)\cdots S\left(j-1\right)$ 在序列S1中至少只出现两次。本文下述前缀标识符，均为对环形序列查找。

3.2. 信息熵理论

3.2.1. 信息熵的定义

信息熵是对信息的一种度量方法，是生物信息学研究中的一种重要工具 [6] [7] [8]。在信息理论中，单符号的离散信息源是最简单的离散信源，一个信源所有可能发送的消息符号构成一个样本空间，用 $X=\left\{x_k|k=1,2,\cdots ,n\right\}$ 表示信源X的样本集合，其中表示 $x_k$ 信源所有发出的消息符号，用 $P_X=\left\{p_k|k=1,2,\cdots ,n\right\}$ 表示概率集合，则 $p_k$ 表示信源发出符号 $x_k$ 的概率。信源X的概率空间表示为

$\left[\begin{array}{c}X\\ P_X\end{array}\right]=\left[\begin{array}{cccc}\begin{array}{l}{x}_1\\ p_1\end{array}& \begin{array}{l}{x}_2\\ p_2\end{array}& \begin{array}{l}\cdots \\ \cdots \end{array}& \begin{array}{l}{x}_n\\ p_n\end{array}\end{array}\right],0\le p_i\le 1,{\displaystyle \underset{i=1}{\overset{n}{\sum }}{p}_i=1}$ (1)

信息熵是信源发出任何一个消息状态所具有的平均信息量 [9]。对于离散无记忆的信源，它的熵定义为 [9]

$H\left(X\right)=-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{p}_i\log p_i}$ (2)

3.2.2. 信息熵熵的性质渐化性

信息熵满足如下条件，该性质称为信息熵的渐化性 [10]。

$H\left(X\right)=H\left(p_1,p_2,\cdots ,p_n\right)=H\left(p_1+p_2,p_3,\cdots ,p_n\right)+\left(p_1+p_2\right)H\left(\frac{p_1}{p_1+p_2},\frac{p_2}{p_1+p_2}\right)$ (3)

$0\le p\left(x_i\right)\le 1,k=1,2,\cdots ,n,{\displaystyle \underset{i=1}{\overset{n}{\sum }}p\left(x_n\right)}=1,p\left(x_1\right)+p\left(x_2\right)>0$

根据上面的等式，发现 $\left(p\left(x_1\right)+p\left(x_2\right)\right)H\left(\frac{p\left(x_1\right)}{p\left(x_1\right)+p\left(x_2\right)},\frac{p\left(x_2\right)}{p\left(x_1\right)+p\left(x_2\right)}\right)$ 是非负数，因此等式右边的

第一项是小于等式左边，则意味着概率合并相加会导致熵变小。也就是说概率合并相加会导致概率分布越不均匀，从而熵变小。因此渐化性反映出概率分布越均匀，熵值越大。

3.3. 相似性度量

任意两条序列之间存在着差异，这些差异可以看作是在进化过程中所产生的突变。本文将两条序列的共同前缀标识符作为它们的共同特征，该标识符在两序列中的位置不同，由此可产生该标识符的位置差，进而有位置差的绝对值。本文认为位置差绝对值的情况越多，意味着进化过程中发生的突变越大，则两序列之间的差异性越大。也就是说位置差绝对值的概率分布越均匀，则两序列的距离越大。于是把该位置差的绝对值看成信源发出的消息符号，每一个绝对值都有相应的出现频率，将该频率看成概率，这些位置差绝对值的信息熵作为两序列之间的距离。

4. 模型

一、提取两序列的共同特征，即找到共同前缀标识符。

(1) 对于任意两条序列S1、S2，分别查找所有位置处的前缀标识符，得到前缀标识符集合，分别记为 $I\left(S1\right)$ 、 $I\left(S2\right)$ ；

(2) 根据前缀标识符集合 $I\left(S1\right)$ 、 $I\left(S2\right)$，得到序列S1、S2共同前缀标识符，记为 $e\left(i\right)$，并记录共同前缀标识符的个数，记为m。其中，

$m=\left|I\left(S1\right)\cap I\left(S2\right)\right|,e\left(i\right)\in I\left(S1\right)\cap I\left(S2\right),i=1,2,\cdots ,m$ (4)

二、计算共同前缀标识符在两序列中位置差的绝对值。

(1) 查找共同前缀标识符 $e\left(i\right)$ 在S1和S2中位置，分别记为 $po{s}^1\left(i\right)$ 、 $po{s}^2\left(i\right)$，其中 $i=1,2,\cdots ,m$ ；

(2) 计算出共同前缀标识符 $e\left(i\right)$ 对应的位置差绝对值，记为 $pos\_dis\left(i\right)$，

$pos\_dis\left(i\right)=\left|po{s}^1\left(i\right)-po{s}^2\left(i\right)\right|,i=1,2,\cdots ,m$ (5)

则所有共同前缀标识符对应的位置差绝对值集合可表达为

$POS\_DIS=\{pos\_dis\left(i\right)|\left|pos\_dis\left(i\right)=po{s}^1\left(i\right)-po{s}^2\left(i\right)\right|,i=1,2,\cdots ,m\}$ (6)

对于任意两个共同前缀标识符 $e\left(i\right)$， $e\left(j\right)$，它们在序列位置差绝对值 $pos\_dis\left(i\right)$ 、 $pos\_dis\left(j\right)$ 可能相等，因此 $POS\_DIS$ 集合为多重集，其中 $i,j\in 1,2,\cdots ,m$。

三、将上述位置差绝对值看成信源发出的符号，得到信源概率空间，计算信息熵，得到两序列间距离。

(1) 根据位置差绝对值对应的多重集 $POS\_DIS$，找到位置差绝对值出现的所有情况，分别记为 $u\_pos\_dis\left(1\right),\cdots ,u\_pos\_dis\left(h\right)$，则共有h中情况；

(2) 统计各类位置差绝对值出现的次数，记为 $n\left(1\right),n\left(2\right),\cdots ,n\left(t\right)$ 则

$n\left(1\right)+\cdots +n\left(h\right)=m$ (7)

(3) 计算各类位置差绝对值出现的频率，记为 $f\left(t\right)$，则

$f\left(t\right)=\frac{n\left(t\right)}{m},t=1,2,\cdots ,h$ (8)

(4) 信源的概率空间为

$\left[\begin{array}{c}X\\ P_X\end{array}\right]=\left[\begin{array}{cccc}\begin{array}{c}u\_pos\_dis\left(1\right)\\ f\left(1\right)\end{array}& \begin{array}{c}u\_pos\_dis\left(2\right)\\ f\left(2\right)\end{array}& \begin{array}{l}\cdots \\ \cdots \end{array}& \begin{array}{c}u\_pos\_dis\left(h\right)\\ f\left(h\right)\end{array}\end{array}\right]$ (9)

(5) 两序列的距离定义为

$dis\left(S1,S2\right)={\displaystyle \underset{t=1}{\overset{h}{\sum }}-f\left(t\right){\log }_{2}f\left(t\right)}$ (10)

信息熵公式中对数的底可以有多种 [10]，本文中计算两序列相似性时取以2为底的对数。

Figure 1. A phylogenetic tree of 70 mammalian mitochondrial DNA sequences by Methods of this article

图1. 本文方法对70条哺乳动物线粒体DNA序列生成的系统发育树

在图1进化树中，蓝色字体表示哺乳动物的超目，黑色字体表示哺乳动物的各个目，进化树的分支代表科，同一目下不同颜色的分支表示不同的科。

Figure 2. A phylogenetic tree of 70 mammalian mitochondrial genomes by The DFT distance of DNA sequences with the 2D numerical mapping [5]

图2. 采用傅里叶分析方法对70条哺乳动物线粒体基因组生成的系统发育树 [5]

5. 结果讨论与分析

本文运用Matlab软件编写代码并实现模型，运用邻接法(Neighbor-joining)构建进化树，运行时间为90秒。

根据生物学分类，哺乳动物分为原生哺乳动物Prototheria、有袋哺乳动物Marsupialia和胎盘哺乳动物Placentalia。其中胎盘哺乳动物又分为灵长总目Euarchontoglires和劳亚兽总目Laurasiatheria。本文的进化树中将原生哺乳动物和有袋哺乳动物分别聚类到各自类别中，但是胎盘哺乳动物聚类欠佳。这是因为在图1进化树中，灵长总目分为四个分支，劳亚兽总目分成两个分支。在图中2进化树中，原生哺乳动物、有袋哺乳动物和胎盘哺乳动物均没有聚类清楚。

从超目的层次分析，本文的进化树将单孔目Monotremata、有袋目Marsupials、非洲兽总目Afrotheria、贫齿目Xenarthra中动物分别聚在各自的超目中，劳亚兽总目Laurasiatheria中只有食虫目中两条序列没有正确聚类，分在劳亚兽总目之外，其余则聚类良好。灵长总目Euarchontoglires由灵长目、啮齿目和兔形目组成，本文中该总目动物聚类欠佳。在图2进化树中，非洲兽总目、贫齿目、劳亚兽总目和灵长总目均没有将序列聚类到各自的类别当中。其中非洲兽总目中序列African elephant序列和aardvark序列与劳亚兽总目中的序列聚在一个节点，劳亚兽总目中hippopotamus序列与灵长总目中序列聚在一个节点，在贫齿目中sloth序列与灵长总目中序列聚在一个节点。

从目的层次分析，本文进化树啮齿目和贫齿目Xenarthra中的皮毛目Pilosa两个目的动物未聚到一起，其余各个目中动物聚类良好。其中鲸偶蹄目有鲸目Cetacea、反刍亚目Ruminantia、猪次目Suina和河马科Hippopotamidae构成，本文这四个目分别聚类良好。鲸偶蹄目在进化关系上也表现良好，与文献 [11] 一致。

在科的层次上分析，本文进化树灵长目中人科中分进来一个长臂猿科动物，鲸偶蹄目中牛科动物分开，其余各个科分别聚类良好。

综上所述，本文提出的序列相似性度量在哺乳动物数据上具有有效性。本文的进化树结果要优于傅里叶分析的进化树结果 [5]。

6. 总结

本文以生物DNA序列以及它们的共同前缀标识符为研究对象，对共同前缀标识符在生物序列中的位置信息进行深入的分析研究，以信息熵作为研究手段，提出一个序列相似性度量方法，建立了一个高效的DNA序列比对模型。依据提出的模型，本文利用70条哺乳动物线粒体DNA序列构建生物进化树，与文献 [5] 提出的DFT方法得到的进化树进行比较，本文的方法则更有优势。

本文提出的方法优点有：(1) 算法的时间复杂度低。(2) 思路简单，易于理解。本文的方法也存在一些不足之处，需要改进完善。该方法只对两序列共同前缀标识符的位置信息进行研究，并没有对共同前缀标识符的长度、以及共同前缀标识符的个数信息研究。

附录

Table 1. Mitochondria DNA sequence information of 70 mammalian genomes

表1. 70条哺乳动物线粒体DNA序列信息