基于流形結(jié)構(gòu)重建的啟動(dòng)子識(shí)別＊

張友新，王立宏

（煙臺(tái)大學(xué)計(jì)算機(jī)學(xué)院，山東 煙臺(tái) 264005）

1 引言

2001 年，Luscombe 對(duì)生物信息學(xué)給出定義[1]：“生物信息學(xué)是根據(jù)分子（從物理化學(xué)角度）和信息技術(shù)（源自應(yīng)用數(shù)學(xué)、計(jì)算機(jī)科學(xué)和統(tǒng)計(jì)學(xué)的原則）的應(yīng)用來(lái)理解和組織與這些分子相關(guān)的大規(guī)模的信息，即生物信息學(xué)是分子生物學(xué)的信息管理系統(tǒng)和諸多實(shí)踐的應(yīng)用”。生物信息學(xué)的主要研究領(lǐng)域涉及基因組學(xué)、蛋白質(zhì)組學(xué)、生物化學(xué)、數(shù)據(jù)挖掘、分子進(jìn)化、分子建模以及算法等[2]，其中DNA 和蛋白質(zhì)是其主要的研究對(duì)象。啟動(dòng)子（Promoter）是DNA 上控制和調(diào)控基因轉(zhuǎn)錄的重要功能元件，因此如何從DNA 序列中快速、準(zhǔn)確識(shí)別出啟動(dòng)子具有重要的研究?jī)r(jià)值。

現(xiàn)有的啟動(dòng)子識(shí)別方法大致可分為三類[3]：（1）基于信號(hào)的方法。其基本思想是通過(guò)識(shí)別不同轉(zhuǎn)錄因子結(jié)合點(diǎn)（TATA 盒、CAAT 盒等元件）的信號(hào)來(lái)識(shí)別啟動(dòng)子，代表算法如Eponine 等[4]。（2）基于內(nèi)容的方法。其基本思想是統(tǒng)計(jì)DNA 中連續(xù)k個(gè)堿基詞分布規(guī)律來(lái)判斷是否為啟動(dòng)子，代表算法如Promoter-Inspector等[5]。（3）基于CpG島的方法。其基本思想是利用在半數(shù)左右的哺乳動(dòng)物啟動(dòng)子附近出現(xiàn)的一段堿基C、G 高密度區(qū)域來(lái)識(shí)別啟動(dòng)子，典型算法如FirstEF等[6]。

基因序列數(shù)據(jù)集通常都具有數(shù)據(jù)量大、維數(shù)高、非線性的特點(diǎn)，因此如何對(duì)基因數(shù)據(jù)進(jìn)行有效的降維處理至關(guān)重要。在我們的前期工作中[7，8]，曾以潛在語(yǔ)義索引的全局模型和潛在語(yǔ)義索引差異模型DLSI（Difference Latent Semantic Indexing）對(duì)基因數(shù)據(jù)進(jìn)行降維處理，其核心技術(shù)是奇異值分解。奇異值分解是線性降維技術(shù)，能有效找出數(shù)據(jù)中的主要關(guān)聯(lián)模式，消除噪聲，因此識(shí)別效果比現(xiàn)有一些啟動(dòng)子識(shí)別算法好[7，8]。Isomap 是非線性降維技術(shù)，降維后能保持樣本間的測(cè)地距離不變，具有較好的流形識(shí)別能力[9]。本文以Isomap為降維工具，研究其流形重建方法對(duì)基因數(shù)據(jù)的適應(yīng)能力，并就參數(shù)調(diào)整和預(yù)處理等幾個(gè)關(guān)鍵問(wèn)題進(jìn)行討論。實(shí)驗(yàn)結(jié)果表明，在設(shè)定的參數(shù)下，對(duì)幾個(gè)大的基因測(cè)試序列的識(shí)別效果均好于線性降維的識(shí)別效果。

2 Isomap及流形結(jié)構(gòu)重建算法

2000年，Tenenbaum[9]在Science上發(fā)表一種非線性降維方法—等距特征映射Isomap（Isometric Mapping），文獻(xiàn)[10]對(duì)其漸進(jìn)收斂性給出了證明。Isomap方法首先在數(shù)據(jù)空間計(jì)算成對(duì)測(cè)地距離，再運(yùn)用多元尺度分析方法[11]把數(shù)據(jù)點(diǎn)從高維輸入空間投影到低維非線性拓?fù)淇臻g中，獲得保持樣本間測(cè)地距離不變的低維流形。Isomap算法具體分為以下三步：

（1）構(gòu)造鄰域圖G：對(duì)于任意兩個(gè)數(shù)據(jù)點(diǎn)i、j，如果j在i的半徑ε 之內(nèi)（或者是i的k 個(gè)最近鄰之一），則連接i 和j，其權(quán)值等于兩點(diǎn)的歐氏距離。

（2）計(jì)算最短路徑：在鄰域圖中計(jì)算任意兩點(diǎn)i和j 之間的最短距離dG（i，j），并以此近似表示拓?fù)淇臻g的測(cè)地距離。

（3）構(gòu)造d 維嵌入：在距離矩陣DG＝｛dG（i，j）｝上，采用經(jīng)典多元尺度分析方法構(gòu)造能保持拓?fù)淇臻g本質(zhì)結(jié)構(gòu)的d 維嵌入空間Y。d 對(duì)應(yīng)殘差下降的拐點(diǎn)，d 之后殘差下降的速度明顯放緩。

該算法自提出以來(lái)得到了研究者的廣泛關(guān)注，由于算法中只涉及到一個(gè)參數(shù)ε（或k），而且算法得到的拓?fù)浣Y(jié)構(gòu)對(duì)該參數(shù)是敏感的[12]，因此人們重點(diǎn)研究有效選擇參數(shù)的方法[13]；為了使結(jié)果的拓?fù)浣Y(jié)構(gòu)更穩(wěn)定，文獻(xiàn)[14]刪除穿越低密度區(qū)域的短路邊，使不同流形之間斷開(kāi)連接。

Isomap方法僅能建立點(diǎn)對(duì)點(diǎn)的降維嵌入，而未能建立高維數(shù)據(jù)流形空間到低維表示空間之間的映射，不能利用訓(xùn)練樣本的降維對(duì)高維空間內(nèi)的測(cè)試樣本進(jìn)行低維表示，從而無(wú)法在低維空間內(nèi)實(shí)現(xiàn)預(yù)測(cè)。文獻(xiàn)[15]分析了Isomap方法的原理，將其本質(zhì)要求概括為連續(xù)性假設(shè)、局部保距假設(shè)和稠密性假設(shè)，并以此為出發(fā)點(diǎn)，推導(dǎo)出了高維流形空間到低維表示空間之間的顯式映射關(guān)系，具體描述如下：

設(shè)高維流形數(shù)據(jù)集為Ωn?Rn，低維表示為Ωd?Rd，映射函數(shù)為g：Ωn→Ωd。訓(xùn)練數(shù)據(jù)集為｛xi｜xi∈Ωn，i＝1，…，m｝，對(duì)應(yīng)的低維表示為｛yi｜yi＝g（xi）∈Ωd，i＝1，…，m｝。測(cè)試數(shù)據(jù)x的低維表示為：

其中，xi是距離x 較近的點(diǎn)，Qi見(jiàn)下面算法描述。

文獻(xiàn)[15]借助映射函數(shù)（1）構(gòu)造出兩種流形結(jié)構(gòu)重建算法：快速算法與穩(wěn)健算法。其中，穩(wěn)健算法考慮到待預(yù)測(cè)樣本的所有ε近鄰信息，通過(guò)映射向量的加權(quán)平均得出待預(yù)測(cè)樣本的映射向量，不易受到噪聲干擾，具有較好的魯棒性。

本文采用的是k最近鄰情況下的穩(wěn)健算法，算法描述如下：

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

上述算法步驟1～步驟3可以看作準(zhǔn)備階段，在存儲(chǔ)空間允許的情況下，可以將計(jì)算結(jié)果保存。步驟4～步驟6是針對(duì)每個(gè)待測(cè)試樣本點(diǎn)進(jìn)行的，步驟6中為保證αj有意義，本文將其分母加1處理。

3 基于流形結(jié)構(gòu)重建的啟動(dòng)子識(shí)別算法

3.1 DNA序列的表示

脫氧核糖核酸（DNA）是主要的遺傳物質(zhì)，核酸的基本結(jié)構(gòu)單位是核苷酸，其組成方式是堿基-戊糖-磷酸。堿基包括嘌呤堿和嘧啶堿兩類，DNA中的嘧啶為胞嘧啶C和腺嘧啶T，嘌呤為腺嘌呤A和鳥(niǎo)嘌呤G。因此，長(zhǎng)鏈形的DNA 序列可以理解為一個(gè)由字母A、C、G、T 組成的字符串。如一段DNA 序列為：

ACTGTAGTCACTGACGTATAATTGACT…

通過(guò)對(duì)序列的組成內(nèi)容進(jìn)行分析，期望找出啟動(dòng)子等功能片斷。啟動(dòng)子是一段核酸序列，它決定DNA 轉(zhuǎn)錄的方向、速度和準(zhǔn)確性。因此，對(duì)啟動(dòng)子的識(shí)別能確定基因轉(zhuǎn)錄起始位點(diǎn)TSS（Transcription Start Site）的大體位置，從而為基因的位置預(yù)測(cè)提供參考。

本文采用詞向量來(lái)表示DNA 序列，這里的“詞”由字符組表示。如固定詞的大小為5個(gè)字符，則可能出現(xiàn)在DNA 序列中的詞共有45＝1 024個(gè)，這些詞形成1 024維的詞向量空間。用長(zhǎng)度為5個(gè)字符的滑動(dòng)窗口截取一個(gè)DNA 序列，每截取一次對(duì)應(yīng)窗口中的特征詞的頻率加1，向下滑動(dòng)窗口步長(zhǎng)為1個(gè)字符，依次截取，這樣就將該序列映射到1 024維的詞向量空間中。

3.2 內(nèi)含子數(shù)據(jù)集的抽樣

除了啟動(dòng)子外，DNA 序列中還含有外顯子（Exon）、內(nèi)含子（Intron）等功能片段，本文的主要目標(biāo)是將一段DNA 序列中的啟動(dòng)子識(shí)別出來(lái)，即與內(nèi)含子、外顯子進(jìn)行區(qū)分。不難發(fā)現(xiàn)，這是一個(gè)三分類問(wèn)題，借助已有的啟動(dòng)子、內(nèi)含子、外顯子訓(xùn)練數(shù)據(jù)集（詳見(jiàn)第4 節(jié)），得出分類模型，對(duì)待測(cè)DNA 序列進(jìn)行測(cè)試，識(shí)別出啟動(dòng)子。

本文采用的分類器是支持向量機(jī)SVM，支持向量機(jī)在各類樣本數(shù)量較為均衡時(shí)分類精度高，因此需要對(duì)這三類數(shù)據(jù)之間進(jìn)行均衡。從樣本數(shù)量上看（詳見(jiàn)第4節(jié)），內(nèi)含子數(shù)量遠(yuǎn)遠(yuǎn)多于啟動(dòng)子和外顯子，因此需要對(duì)內(nèi)含子集合進(jìn)行抽樣，以保留有代表性的樣本?；诮弬鞑サ木垲愃惴ˋP（Affinity Propagation）是2007 年Frey B J提出的[16]，AP能得出有意義的代表點(diǎn)集合，因而可以作為有效抽樣的工具使用。關(guān)于AP 算法的實(shí)現(xiàn)細(xì)節(jié)見(jiàn)文獻(xiàn)[16]，偏向參數(shù)通常設(shè)置為相似度的中位數(shù)。為了獲取和啟動(dòng)子、外顯子集合大小相當(dāng)?shù)膬?nèi)含子集合，本文首先將內(nèi)含子數(shù)據(jù)經(jīng)3.4節(jié)的詞向量歸一化后，所有內(nèi)含子之間的相似度按從小到大排序，設(shè)置偏向參數(shù)為相似度序列0.96位置的值，經(jīng)過(guò)AP抽樣得到內(nèi)含子集合的大小為769。

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

3.3 訓(xùn)練數(shù)據(jù)的Isomap降維

為了觀察數(shù)據(jù)的分布情況，本文將這三類數(shù)據(jù)混合，共計(jì)565＋769＋890＝2 224 個(gè)序列，1 024維，采用Isomap 降維，得到各類數(shù)據(jù)三維分布圖如圖2a所示，得到的殘差圖如圖2b所示。

Figure 2 Result of isomap dimensionality reduction圖2 Isomap降維結(jié)果

從圖2a中可以看出，三類數(shù)據(jù)之間并不是線性可分的，數(shù)據(jù)之間混疊情況較嚴(yán)重。長(zhǎng)DNA 序列中識(shí)別啟動(dòng)子即以此作為訓(xùn)練樣本，找出測(cè)試數(shù)據(jù)中的啟動(dòng)子。從圖2b 可以看出，從維數(shù)d＝5開(kāi)始，殘差下降速度放緩，本文取d＝6。

3.4 重建算法的準(zhǔn)備條件

公式（1）推導(dǎo)過(guò)程中采用近似的方法，對(duì)O（ε2）的項(xiàng)忽略不計(jì)，因此在使用過(guò)程中，為了保證近似的效果，需要盡量減小ε的值。本文考慮的是k鄰域情況下，此時(shí)應(yīng)減小鄰域內(nèi)各點(diǎn)之間的距離。前期實(shí)驗(yàn)表明，如果不采取任何措施，直接使用流形重建算法對(duì)啟動(dòng)子的識(shí)別效果是非常差的。這里采用的方法有：特征選擇和詞向量歸一化。通過(guò)特征選擇減少維數(shù)，使兩點(diǎn)之間歐氏距離減小；通過(guò)向量歸一化減小每個(gè)數(shù)據(jù)值，也能減小歐氏距離。

3.4.1 特征選擇

由于這三類數(shù)據(jù)的詞向量空間完全相同，而各類數(shù)據(jù)在詞向量空間的分布不同，因此本文借助詞頻統(tǒng)計(jì)與過(guò)濾等方法，找出每類數(shù)據(jù)相關(guān)程度較高的詞向量子空間，以此完成初步特征選擇。

以啟動(dòng)子為例，詞頻統(tǒng)計(jì)與過(guò)濾方法如下：

（1）統(tǒng)計(jì)詞頻，生成特征詞-序列矩陣。用3.1節(jié)所述的滑動(dòng)窗口方法將每個(gè)啟動(dòng)子序列映射到1 024維的詞向量空間中，565 個(gè)啟動(dòng)子對(duì)應(yīng)矩陣XP的大小為1 024×565。

類似地得到內(nèi)含子和外顯子對(duì)應(yīng)的矩陣XI和XE，對(duì)應(yīng)的特征詞集合WI和WE。為了進(jìn)一步區(qū)分啟動(dòng)子和內(nèi)含子，定義啟動(dòng)子和內(nèi)含子的差異特征詞集合WPI為WP和WI的異或，記錄那些只在WP或只在WI中出現(xiàn)的特征詞。同樣，定義啟動(dòng)子和外顯子的差異特征詞集合WPE為WP和WE的異或。

3.4.2 詞向量歸一化

3.5 基于流形重建的啟動(dòng)子識(shí)別算法

本文采用支持向量機(jī)LIBSVM （http：／／www.csie.ntu.edu.tw／～cjlin／）進(jìn)行分類，分別建立啟動(dòng)子-內(nèi)含子和啟動(dòng)子-外顯子分類器，當(dāng)且僅當(dāng)這兩個(gè)分類器都認(rèn)定是啟動(dòng)子時(shí)算法才判斷為啟動(dòng)子，如圖3所示。

具體流程如下：

（1）采用3.1節(jié)所述方法將啟動(dòng)子、內(nèi)含子、外顯子訓(xùn)練序列和待測(cè)試序列（截取長(zhǎng)300字符的小段，滑動(dòng)5個(gè)字符，得到下一小段）都映射到1 024維的詞向量空間中，得到矩陣XP、XI、XE和Test。記錄差異特征詞集合WPI和WPE。

Figure 3 Classification method圖3 分類方法

（2）將合并矩陣[XP，XI]與WPI所對(duì)應(yīng)的詞項(xiàng)形成的子矩陣作為訓(xùn)練數(shù)據(jù)，Test與WPI所對(duì)應(yīng)的詞項(xiàng)形成的子矩陣作為測(cè)試數(shù)據(jù)，輸入穩(wěn)健的流形結(jié)構(gòu)映射算法進(jìn)行降維（見(jiàn)第2節(jié)），輸出[XP，XI]的低維嵌入YPI和Test 的低維嵌入Yout1。矩陣[XP，XE]的處理方法相同，輸出為YPE和Yout2。

（3）將YPI附加上類別信息后輸入支持向量機(jī)工具LIBSVM 中進(jìn)行訓(xùn)練，得到分類模型YPI.model；同樣地，YPE附加上類別信息后進(jìn)行訓(xùn)練，得到分類模型YPE.model。

（4）在LIBSVM 中，用YPI.model對(duì)Yout1進(jìn)行預(yù)測(cè)分類，用YPE.model對(duì)Yout2進(jìn)行預(yù)測(cè)分類。

4 實(shí)驗(yàn)結(jié)果與分析

基于流形結(jié)構(gòu)重建的啟動(dòng)子識(shí)別算法分為訓(xùn)練和測(cè)試兩個(gè)階段，訓(xùn)練數(shù)據(jù)集為美國(guó)勞倫斯伯克利國(guó)家實(shí)驗(yàn)室和加州大學(xué)圣克魯茲分校提供的一組用于基因識(shí)別算法的公共數(shù)據(jù)集，下載地址為http：／／www.fruitfly.org／seq＿tools／datasets／Human／，數(shù)據(jù)集基本情況如表1所示。測(cè)試數(shù)據(jù)采用與文獻(xiàn)[3]和文獻(xiàn)[7]相同的六個(gè)DNA 測(cè)試序列，各序列的基本情況如表2所示。

Table 1 Basic information of training datasets表1 訓(xùn)練數(shù)據(jù)集基本情況

Table 2 Basic information of test datasets表2 測(cè)試序列基本情況

4.1 參數(shù)設(shè)定

Figure 4 Relationship between and the number of feature words圖4 與特征詞項(xiàng)數(shù)目的關(guān)系

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

對(duì)于啟動(dòng)子識(shí)別結(jié)果的處理，采用與文獻(xiàn)[7]相同的方法，即在某一位置附近連續(xù)出現(xiàn)的識(shí)別結(jié)果視為同一個(gè)啟動(dòng)子，并按認(rèn)定次數(shù)降序排列，在轉(zhuǎn)錄起始點(diǎn)的[-2 000，500]的認(rèn)定結(jié)果有效。采用準(zhǔn)確率和查全率來(lái)比較和評(píng)價(jià)基于流形結(jié)構(gòu)重建的啟動(dòng)子識(shí)別算法的識(shí)別效果，如表3所示。除本文算法外，其余實(shí)驗(yàn)結(jié)果來(lái)自文獻(xiàn)[7]。DLSI是采用線性降維技術(shù)得出的識(shí)別結(jié)果。

從表3 可以看出，對(duì)于測(cè)試序列AF017257、AC002368和D87675，本文算法的結(jié)果與其它算法的最好效果相同，達(dá)到100%的查全率和正確率；對(duì)于測(cè)試序列AF146793、L44140，本文算法的結(jié)果與其它算法的最好效果正確識(shí)別的啟動(dòng)子個(gè)數(shù)相同，而錯(cuò)誤識(shí)別的啟動(dòng)子個(gè)數(shù)下降；對(duì)于測(cè)試序列AC002397，本文算法得到的正確啟動(dòng)子更多，查全率最高，而錯(cuò)誤識(shí)別的啟動(dòng)子并未明顯增加。因此，綜合各項(xiàng)指標(biāo)，本文算法與文獻(xiàn)[7]相比，實(shí)驗(yàn)結(jié)果較好。

5 結(jié)束語(yǔ)

本文以文本分析的方法研究了生物信息學(xué)中的啟動(dòng)子識(shí)別問(wèn)題。由于DNA 序列表示為詞向量空間后得到的是高維數(shù)據(jù)，對(duì)數(shù)據(jù)的有效降維是不可避免的問(wèn)題。本文采用Isomap 方法為其降維，但I(xiàn)somap方法只能降維不能預(yù)測(cè)。為了能預(yù)測(cè)長(zhǎng)DNA 序列中啟動(dòng)子的位置，本文借鑒流形結(jié)構(gòu)重建算法對(duì)訓(xùn)練數(shù)據(jù)降維后的信息加以利用，從而完成預(yù)測(cè)。直接使用重建算法不能滿足其稠密性假設(shè)，故采用特征選擇和向量歸一化等方法減小鄰域內(nèi)點(diǎn)之間的距離，使預(yù)測(cè)結(jié)果可信。通過(guò)對(duì)六個(gè)長(zhǎng)DNA 序列進(jìn)行測(cè)試，本文實(shí)驗(yàn)方法好于文獻(xiàn)結(jié)果。

Table 3 Comparison of experiment results表3 實(shí)驗(yàn)結(jié)果對(duì)比

[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.

附中文參考文獻(xiàn)：

[1]鐘楊，張亮，趙瓊.簡(jiǎn)明生物信息學(xué)[M].北京：高等教育出版社，2001.

[8]秦洋，王立宏，武栓虎，等.啟動(dòng)子的潛在語(yǔ)義索引差異識(shí)別算法[J].煙臺(tái)大學(xué)學(xué)報(bào)：自然科學(xué)與工程版，2010，2.（3）：211-216.

[14]邵超，黃厚寬，趙連偉.一種更具拓?fù)浞€(wěn)定性的isomap算法[J].軟件學(xué)報(bào)，2007，18（4）：869-877.

[15]孟德宇，徐晨，徐宗本.基于Isomap的流形結(jié)構(gòu)重建方法[J].計(jì)算機(jī)學(xué)報(bào)，2010，33（3）：545-555.