張友新，王立宏

（煙臺大學計算機學院，山東 煙臺 264005）

1 引言

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

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

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

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

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

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

（2）計算最短路徑：在鄰域圖中計算任意兩點i和j 之間的最短距離dG（i，j），并以此近似表示拓撲空間的測地距離。

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

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

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

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

其中，xi是距離x 較近的點，Qi見下面算法描述。

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

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

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

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

3 基于流形結(jié)構重建的啟動子識別算法

3.1 DNA序列的表示

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

ACTGTAGTCACTGACGTATAATTGACT…

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

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

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

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

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

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

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

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

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

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

3.4 重建算法的準備條件

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

3.4.1 特征選擇

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

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

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

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

3.4.2 詞向量歸一化

3.5 基于流形重建的啟動子識別算法

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

具體流程如下：

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

Figure 3 Classification method圖3 分類方法

（2）將合并矩陣[XP，XI]與WPI所對應的詞項形成的子矩陣作為訓練數(shù)據(jù)，Test與WPI所對應的詞項形成的子矩陣作為測試數(shù)據(jù)，輸入穩(wěn)健的流形結(jié)構映射算法進行降維（見第2節(jié)），輸出[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 實驗結(jié)果與分析

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

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

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

4.1 參數(shù)設定

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

4.2 實驗結(jié)果

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

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

5 結(jié)束語

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

Table 3 Comparison of experiment results表3 實驗結(jié)果對比

[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]邵超，黃厚寬，趙連偉.一種更具拓撲穩(wěn)定性的isomap算法[J].軟件學報，2007，18（4）：869-877.

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