張友新,王立宏
(煙臺(tái)大學(xué)計(jì)算機(jī)學(xué)院,山東 煙臺(tái) 264005)
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í)別效果。
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處理。
脫氧核糖核酸(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維的詞向量空間中。
除了啟動(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 排序后的相似度序列
為了觀察數(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。
公式(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 詞向量歸一化
本文采用支持向量機(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è)分類。
基于流形結(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è)試序列基本情況
Figure 4 Relationship between and the number of feature words圖4 與特征詞項(xiàng)數(shù)目的關(guān)系
對(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é)果較好。
本文以文本分析的方法研究了生物信息學(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.