基于流形结构重建的启动子识别＊

张友新，王立宏

（烟台大学计算机学院，山东 烟台 264005）

1 引言

2001 年，Luscombe 对生物信息学给出定义[1]：“生物信息学是根据分子（从物理化学角度）和信息技术（源自应用数学、计算机科学和统计学的原则）的应用来理解和组织与这些分子相关的大规模的信息，即生物信息学是分子生物学的信息管理系统和诸多实践的应用”。生物信息学的主要研究领域涉及基因组学、蛋白质组学、生物化学、数据挖掘、分子进化、分子建模以及算法等[2]，其中DNA 和蛋白质是其主要的研究对象。启动子（Promoter）是DNA 上控制和调控基因转录的重要功能元件，因此如何从DNA 序列中快速、准确识别出启动子具有重要的研究价值。

现有的启动子识别方法大致可分为三类[3]：（1）基于信号的方法。其基本思想是通过识别不同转录因子结合点（TATA 盒、CAAT 盒等元件）的信号来识别启动子，代表算法如Eponine 等[4]。（2）基于内容的方法。其基本思想是统计DNA 中连续k个碱基词分布规律来判断是否为启动子，代表算法如Promoter-Inspector等[5]。（3）基于CpG岛的方法。其基本思想是利用在半数左右的哺乳动物启动子附近出现的一段碱基C、G 高密度区域来识别启动子，典型算法如FirstEF等[6]。

基因序列数据集通常都具有数据量大、维数高、非线性的特点，因此如何对基因数据进行有效的降维处理至关重要。在我们的前期工作中[7，8]，曾以潜在语义索引的全局模型和潜在语义索引差异模型DLSI（Difference Latent Semantic Indexing）对基因数据进行降维处理，其核心技术是奇异值分解。奇异值分解是线性降维技术，能有效找出数据中的主要关联模式，消除噪声，因此识别效果比现有一些启动子识别算法好[7，8]。Isomap 是非线性降维技术，降维后能保持样本间的测地距离不变，具有较好的流形识别能力[9]。本文以Isomap为降维工具，研究其流形重建方法对基因数据的适应能力，并就参数调整和预处理等几个关键问题进行讨论。实验结果表明，在设定的参数下，对几个大的基因测试序列的识别效果均好于线性降维的识别效果。

2 Isomap及流形结构重建算法

2000年，Tenenbaum[9]在Science上发表一种非线性降维方法—等距特征映射Isomap（Isometric Mapping），文献[10]对其渐进收敛性给出了证明。Isomap方法首先在数据空间计算成对测地距离，再运用多元尺度分析方法[11]把数据点从高维输入空间投影到低维非线性拓扑空间中，获得保持样本间测地距离不变的低维流形。Isomap算法具体分为以下三步：

（1）构造邻域图G：对于任意两个数据点i、j，如果j在i的半径ε 之内（或者是i的k 个最近邻之一），则连接i 和j，其权值等于两点的欧氏距离。

（2）计算最短路径：在邻域图中计算任意两点i和j 之间的最短距离dG（i，j），并以此近似表示拓扑空间的测地距离。

（3）构造d 维嵌入：在距离矩阵DG＝｛dG（i，j）｝上，采用经典多元尺度分析方法构造能保持拓扑空间本质结构的d 维嵌入空间Y。d 对应残差下降的拐点，d 之后残差下降的速度明显放缓。

该算法自提出以来得到了研究者的广泛关注，由于算法中只涉及到一个参数ε（或k），而且算法得到的拓扑结构对该参数是敏感的[12]，因此人们重点研究有效选择参数的方法[13]；为了使结果的拓扑结构更稳定，文献[14]删除穿越低密度区域的短路边，使不同流形之间断开连接。

Isomap方法仅能建立点对点的降维嵌入，而未能建立高维数据流形空间到低维表示空间之间的映射，不能利用训练样本的降维对高维空间内的测试样本进行低维表示，从而无法在低维空间内实现预测。文献[15]分析了Isomap方法的原理，将其本质要求概括为连续性假设、局部保距假设和稠密性假设，并以此为出发点，推导出了高维流形空间到低维表示空间之间的显式映射关系，具体描述如下：

设高维流形数据集为Ωn⊂Rn，低维表示为Ωd⊂Rd，映射函数为g：Ωn→Ωd。训练数据集为｛xi｜xi∈Ωn，i＝1，…，m｝，对应的低维表示为｛yi｜yi＝g（xi）∈Ωd，i＝1，…，m｝。测试数据x的低维表示为：

其中，xi是距离x 较近的点，Qi见下面算法描述。

文献[15]借助映射函数（1）构造出两种流形结构重建算法：快速算法与稳健算法。其中，稳健算法考虑到待预测样本的所有ε近邻信息，通过映射向量的加权平均得出待预测样本的映射向量，不易受到噪声干扰，具有较好的鲁棒性。

本文采用的是k最近邻情况下的稳健算法，算法描述如下：

输出：任意x∈Ωn的映射向量y ∈Ωd。

上述算法步骤1～步骤3可以看作准备阶段，在存储空间允许的情况下，可以将计算结果保存。步骤4～步骤6是针对每个待测试样本点进行的，步骤6中为保证αj有意义，本文将其分母加1处理。

3 基于流形结构重建的启动子识别算法

3.1 DNA序列的表示

脱氧核糖核酸（DNA）是主要的遗传物质，核酸的基本结构单位是核苷酸，其组成方式是碱基-戊糖-磷酸。碱基包括嘌呤碱和嘧啶碱两类，DNA中的嘧啶为胞嘧啶C和腺嘧啶T，嘌呤为腺嘌呤A和鸟嘌呤G。因此，长链形的DNA 序列可以理解为一个由字母A、C、G、T 组成的字符串。如一段DNA 序列为：

ACTGTAGTCACTGACGTATAATTGACT…

通过对序列的组成内容进行分析，期望找出启动子等功能片断。启动子是一段核酸序列，它决定DNA 转录的方向、速度和准确性。因此，对启动子的识别能确定基因转录起始位点TSS（Transcription Start Site）的大体位置，从而为基因的位置预测提供参考。

本文采用词向量来表示DNA 序列，这里的“词”由字符组表示。如固定词的大小为5个字符，则可能出现在DNA 序列中的词共有45＝1 024个，这些词形成1 024维的词向量空间。用长度为5个字符的滑动窗口截取一个DNA 序列，每截取一次对应窗口中的特征词的频率加1，向下滑动窗口步长为1个字符，依次截取，这样就将该序列映射到1 024维的词向量空间中。

3.2 内含子数据集的抽样

除了启动子外，DNA 序列中还含有外显子（Exon）、内含子（Intron）等功能片段，本文的主要目标是将一段DNA 序列中的启动子识别出来，即与内含子、外显子进行区分。不难发现，这是一个三分类问题，借助已有的启动子、内含子、外显子训练数据集（详见第4 节），得出分类模型，对待测DNA 序列进行测试，识别出启动子。

本文采用的分类器是支持向量机SVM，支持向量机在各类样本数量较为均衡时分类精度高，因此需要对这三类数据之间进行均衡。从样本数量上看（详见第4节），内含子数量远远多于启动子和外显子，因此需要对内含子集合进行抽样，以保留有代表性的样本。基于近邻传播的聚类算法AP（Affinity Propagation）是2007 年Frey B J提出的[16]，AP能得出有意义的代表点集合，因而可以作为有效抽样的工具使用。关于AP 算法的实现细节见文献[16]，偏向参数通常设置为相似度的中位数。为了获取和启动子、外显子集合大小相当的内含子集合，本文首先将内含子数据经3.4节的词向量归一化后，所有内含子之间的相似度按从小到大排序，设置偏向参数为相似度序列0.96位置的值，经过AP抽样得到内含子集合的大小为769。

Figure 1 Sorted similarity sequences图1 排序后的相似度序列

3.3 训练数据的Isomap降维

为了观察数据的分布情况，本文将这三类数据混合，共计565＋769＋890＝2 224 个序列，1 024维，采用Isomap 降维，得到各类数据三维分布图如图2a所示，得到的残差图如图2b所示。

Figure 2 Result of isomap dimensionality reduction图2 Isomap降维结果

从图2a中可以看出，三类数据之间并不是线性可分的，数据之间混叠情况较严重。长DNA 序列中识别启动子即以此作为训练样本，找出测试数据中的启动子。从图2b 可以看出，从维数d＝5开始，残差下降速度放缓，本文取d＝6。

3.4 重建算法的准备条件

公式（1）推导过程中采用近似的方法，对O（ε2）的项忽略不计，因此在使用过程中，为了保证近似的效果，需要尽量减小ε的值。本文考虑的是k邻域情况下，此时应减小邻域内各点之间的距离。前期实验表明，如果不采取任何措施，直接使用流形重建算法对启动子的识别效果是非常差的。这里采用的方法有：特征选择和词向量归一化。通过特征选择减少维数，使两点之间欧氏距离减小；通过向量归一化减小每个数据值，也能减小欧氏距离。

3.4.1 特征选择

由于这三类数据的词向量空间完全相同，而各类数据在词向量空间的分布不同，因此本文借助词频统计与过滤等方法，找出每类数据相关程度较高的词向量子空间，以此完成初步特征选择。

以启动子为例，词频统计与过滤方法如下：

（1）统计词频，生成特征词-序列矩阵。用3.1节所述的滑动窗口方法将每个启动子序列映射到1 024维的词向量空间中，565 个启动子对应矩阵XP的大小为1 024×565。

类似地得到内含子和外显子对应的矩阵XI和XE，对应的特征词集合WI和WE。为了进一步区分启动子和内含子，定义启动子和内含子的差异特征词集合WPI为WP和WI的异或，记录那些只在WP或只在WI中出现的特征词。同样，定义启动子和外显子的差异特征词集合WPE为WP和WE的异或。

3.4.2 词向量归一化

3.5 基于流形重建的启动子识别算法

本文采用支持向量机LIBSVM （http：／／www.csie.ntu.edu.tw／～cjlin／）进行分类，分别建立启动子-内含子和启动子-外显子分类器，当且仅当这两个分类器都认定是启动子时算法才判断为启动子，如图3所示。

具体流程如下：

（1）采用3.1节所述方法将启动子、内含子、外显子训练序列和待测试序列（截取长300字符的小段，滑动5个字符，得到下一小段）都映射到1 024维的词向量空间中，得到矩阵XP、XI、XE和Test。记录差异特征词集合WPI和WPE。

Figure 3 Classification method图3 分类方法

（2）将合并矩阵[XP，XI]与WPI所对应的词项形成的子矩阵作为训练数据，Test与WPI所对应的词项形成的子矩阵作为测试数据，输入稳健的流形结构映射算法进行降维（见第2节），输出[XP，XI]的低维嵌入YPI和Test 的低维嵌入Yout1。矩阵[XP，XE]的处理方法相同，输出为YPE和Yout2。

（3）将YPI附加上类别信息后输入支持向量机工具LIBSVM 中进行训练，得到分类模型YPI.model；同样地，YPE附加上类别信息后进行训练，得到分类模型YPE.model。

（4）在LIBSVM 中，用YPI.model对Yout1进行预测分类，用YPE.model对Yout2进行预测分类。

4 实验结果与分析

基于流形结构重建的启动子识别算法分为训练和测试两个阶段，训练数据集为美国劳伦斯伯克利国家实验室和加州大学圣克鲁兹分校提供的一组用于基因识别算法的公共数据集，下载地址为http：／／www.fruitfly.org／seq＿tools／datasets／Human／，数据集基本情况如表1所示。测试数据采用与文献[3]和文献[7]相同的六个DNA 测试序列，各序列的基本情况如表2所示。

Table 1 Basic information of training datasets表1 训练数据集基本情况

Table 2 Basic information of test datasets表2 测试序列基本情况

4.1 参数设定

Figure 4 Relationship between and the number of feature words图4 与特征词项数目的关系

4.2 实验结果

对于启动子识别结果的处理，采用与文献[7]相同的方法，即在某一位置附近连续出现的识别结果视为同一个启动子，并按认定次数降序排列，在转录起始点的[-2 000，500]的认定结果有效。采用准确率和查全率来比较和评价基于流形结构重建的启动子识别算法的识别效果，如表3所示。除本文算法外，其余实验结果来自文献[7]。DLSI是采用线性降维技术得出的识别结果。

从表3 可以看出，对于测试序列AF017257、AC002368和D87675，本文算法的结果与其它算法的最好效果相同，达到100%的查全率和正确率；对于测试序列AF146793、L44140，本文算法的结果与其它算法的最好效果正确识别的启动子个数相同，而错误识别的启动子个数下降；对于测试序列AC002397，本文算法得到的正确启动子更多，查全率最高，而错误识别的启动子并未明显增加。因此，综合各项指标，本文算法与文献[7]相比，实验结果较好。

5 结束语

本文以文本分析的方法研究了生物信息学中的启动子识别问题。由于DNA 序列表示为词向量空间后得到的是高维数据，对数据的有效降维是不可避免的问题。本文采用Isomap 方法为其降维，但Isomap方法只能降维不能预测。为了能预测长DNA 序列中启动子的位置，本文借鉴流形结构重建算法对训练数据降维后的信息加以利用，从而完成预测。直接使用重建算法不能满足其稠密性假设，故采用特征选择和向量归一化等方法减小邻域内点之间的距离，使预测结果可信。通过对六个长DNA 序列进行测试，本文实验方法好于文献结果。

Table 3 Comparison of experiment results表3 实验结果对比

[1]Zhong Liang，Zhang Liang，Zhao Qiong.An introduction to bioinformatics[M].Beijing：Higher Education Press，2001.（in Chinese）

[2]Trifonov E N.Earliest pages of bioinformatics[J].Bioinformatics，2000，16（1）：5-9.

[3]Wu Shuan-hu，Xie Xu-dong，Wee-Chung L A，et al.Eukaryotic promoter prediction based on relative entropy and positional information[J].Physical Review E，2007，75（4）：1-7.

[4]Down T A，Hubbard T J.Computational detection and location of transcription start sites in mammalian genomic DNA[J].Henome Research，2002（12）：458-461.

[5]Scherf M，Klingenhoff A，Werner T.Highly specific localization of promoter regions in large genomic sequences by PromoterInspector：a novel context analysis approach[J].Journal of Molecular Biology，2000，2.7（3）：599-606.

[6]Davuluri R.V，Grosse I，Zhang M Q.Computational identification of promoters and first exons in the human genome[J].Nature Genetics，2001（29）：412-417.

[7]Wang Li-hong，Zhao Xian-jia，Qin Yang.Latent semantic indexing in eukaryotic promoter recognition[J].Journal of Computational Information Systems，2010，6（5）：1509-1516.

[8]Qin Yang，Wang Li-hong，Wu Shuan-hu，et al.Genomic promoter recognition algorithm based on difference latent semantic indexing model[J].Journal of Yantai University：Natural Science and Engineering Edition，2010，2.（3）：211-216.（in Chinese）

[9]Tenenbaum J B，de Silva V，Langford J C.A global geometric framework for nonlinear dimensionality reduction[J].Science，2000，2.0（5500）：2319-2323.

[10]Bernstein M，de Silva V，Langford J C，et al.Graph approximations to geodesics on embedded manifolds[EB／OL].[2000-01-01].Technical Report，Stanford University，2000.http：／／isomap.stanford.edu／BdSLT.pdf.

[11]Cox T，Cox M.Multidimensional scaling[M].London：Chapman ＆ Hall，1994.

[12]Balasubramanian M，Schwartz E L，et al.The Isomap algorithm and topological stability[J].Science，2002，2.5（5552）：7.

[13]Shao C，Huang H K.Selection of the optimal parameter value for the isomap algorithm[C]∥Proc of the 4th Mexican Int'l Conf on Artificial Intelligence，2005：396-404.

[14]Shao Chao，Huang Hou-kuan，Zhao Lian-wei.A more to-pologically stable isomap algorithm[J].Journal of Software，2007，18（4）：869-877.（in Chinese）

[15]Meng De-yu，Xu Chen，Xu Zong-ben.A new manifold reconstruction method based on Isomap[J].Chinese Journal of Computers，2010，33（3），545-555.（in Chinese）

[16]Frey B J，Dueck D.Clustering by passing messages between data points[J].Science，2007，315（5814）：927-976.

附中文参考文献：

[1]钟杨，张亮，赵琼.简明生物信息学[M].北京：高等教育出版社，2001.

[8]秦洋，王立宏，武栓虎，等.启动子的潜在语义索引差异识别算法[J].烟台大学学报：自然科学与工程版，2010，2.（3）：211-216.

[14]邵超，黄厚宽，赵连伟.一种更具拓扑稳定性的isomap算法[J].软件学报，2007，18（4）：869-877.

[15]孟德宇，徐晨，徐宗本.基于Isomap的流形结构重建方法[J].计算机学报，2010，33（3）：545-555.