郭海濤 霍紅衛(wèi) 于強

SegHMC:一種基于Segmental HMM模型的順式調(diào)控模塊識別算法

郭海濤1霍紅衛(wèi)1于強1

順式調(diào)控模塊(Cis-regulatory module,CRM)在真核生物基因的轉(zhuǎn)錄調(diào)控中起著重要作用,識別順式調(diào)控模塊是當(dāng)前計算生物學(xué)的一個重要課題.雖然當(dāng)前有許多計算方法用于識別順式調(diào)控模塊,但識別準(zhǔn)確率仍有待進一步提高.將順式調(diào)控模塊的多種特征信息結(jié)合在一起,有助于提高識別順式調(diào)控模塊的準(zhǔn)確率.基于此,本文提出了一種識別順式調(diào)控模塊的算法SegHMC(Segmental HMM model for discovery of cis-regulatory module).該算法建立了一種關(guān)于順式調(diào)控模塊識別問題的Segmental HMM模型,進一步擴展了順式調(diào)控模塊調(diào)控結(jié)構(gòu)(或調(diào)控語法)的表示,不僅將順式調(diào)控模塊表示為模體(Motif)的組合,還進一步將模體共同出現(xiàn)的頻率、模體順序偏好以及順式調(diào)控模塊中相鄰模體間的距離分布等特征引入到順式調(diào)控模塊的調(diào)控語法中.在模擬數(shù)據(jù)集和真實生物數(shù)據(jù)集上的實驗結(jié)果表明,本文方法識別順式調(diào)控模塊的準(zhǔn)確率顯著優(yōu)于當(dāng)前的主要方法.

基因的轉(zhuǎn)錄調(diào)控,模體,Segmental HMM,順式調(diào)控模塊識別

在基因的表達調(diào)控系統(tǒng)中,轉(zhuǎn)錄因子(Transcription factor,TF)通過與所調(diào)控基因附近被稱為轉(zhuǎn)錄因子結(jié)合位點(Transcription factor binding site,TFBS)或模體(Motif)的特定DNA序列片段相結(jié)合,來啟動基因的轉(zhuǎn)錄調(diào)控[1?2].在真核生物中,多個轉(zhuǎn)錄因子對基因的轉(zhuǎn)錄調(diào)控,并不是孤立進行的,而是轉(zhuǎn)錄因子之間或者各個轉(zhuǎn)錄因子與它們的模體之間通過一系列的時空交互來實施更復(fù)雜、更精確的轉(zhuǎn)錄調(diào)控.在被調(diào)控基因的調(diào)控區(qū)(Transcriptional regulatory region)中,模體非均勻地聚集為一系列的稱為順式調(diào)控模塊(Cis-regulatory module,CRM)的離散區(qū)域,如啟動子(Promoter)、增強子(Enhancer)、沉寂子(Silencer)、絕緣子(Insulator)等.結(jié)合這些順式調(diào)控模塊的轉(zhuǎn)錄因子通過相互協(xié)作、相互競爭,激活或抑制所調(diào)控基因的轉(zhuǎn)錄表達.一個順式調(diào)控模塊包含單個或多個轉(zhuǎn)錄因子的多個模體實例,其長度通常約為幾百到幾千個堿基對(Base pair,bp).一個真核生物基因調(diào)控區(qū)序列可能的調(diào)控結(jié)構(gòu)的簡單例子如圖1所示.

圖1 順式調(diào)控模塊結(jié)構(gòu)示意圖(順式調(diào)控模塊是包含多個轉(zhuǎn)錄因子相應(yīng)模體的序列區(qū);模體的方向、模體間的間隔距離、模體間的相互關(guān)系可能包含了給定順式調(diào)控模塊的重要性質(zhì).)Fig.1 The structure discription of cis-regulatory modules(A cis-regulatory module is a sequence region that contains multiple motifs of multiple transcription factors;motif orientation,the interval distance between motifs and their cooperation relationship may imply the important regulatory properties of the cis-regulatory module.)

識別順式調(diào)控模塊是理解基因轉(zhuǎn)錄調(diào)控分子機制的基礎(chǔ),同時也是構(gòu)建基因調(diào)控網(wǎng)絡(luò)[3?4]的關(guān)鍵步驟.此外,識別具有特定調(diào)控功能的順式調(diào)控模塊對疾病機理的研究也有重要的意義.許多疾病的發(fā)生都與基因的異常表達有關(guān),調(diào)控基因表達的順式調(diào)控模塊發(fā)生變異是造成基因異常表達的主因.有證據(jù)表明特定順式調(diào)控模塊中協(xié)作調(diào)控元素的破壞,可以導(dǎo)致畸形和疾病;例如,Kleinjan等[5]發(fā)現(xiàn)PAX6的任何遠(yuǎn)端調(diào)控元素的缺失都會改變其表達水平,從而造成先天性眼球畸形、無虹膜以及大腦缺陷等疾病.

通過生物實驗,例如高通量檢測與順式調(diào)控模塊相關(guān)的表觀遺傳標(biāo)記特征[6],可以識別順式調(diào)控模塊,但這種方法費時費力代價較大,并且許多時候受限于實驗條件而很難實施.因此,使用計算方法直接從DNA序列中識別順式調(diào)控模塊已成為一個非常有吸引力的手段.

然而,使用計算方法識別順式調(diào)控模塊也存在著許多挑戰(zhàn):1)同一個順式調(diào)控模塊不同實例內(nèi)的模體排列順序并不完全相同,但也并非完全無序的.此外,模塊內(nèi)的模體之間的距離也不確定,即同一模塊不同實例內(nèi)相同的兩相鄰模體間的距離也都不相同.因此,很難確定性地刻畫這種結(jié)構(gòu).2)真核生物的調(diào)控區(qū)通常很長,構(gòu)成順式調(diào)控模塊的模體通常較短且存在退化,根據(jù)已知的模體或者借助于現(xiàn)有的模體庫(如TRANSFAC[7]、JASPAR[8])直接搜索,會找出大量的假陽性匹配,所以很難通過直接搜索相關(guān)模體的方式來識別包含這些模體的順式調(diào)控模塊.

至今,已有多種用于識別順式調(diào)控模塊的模型和方法[6,9?27].為了識別真核生物基因的順式調(diào)控模塊,不同方法利用順式調(diào)控模塊的不同特征(如模體的聚集和物種間的保守性),使用不同搜索策略.

其中一類方法基于窗聚集,利用模體傾向于聚集的特性來搜索順式調(diào)控模塊.這類方法用最簡單的方式表示順式調(diào)控模塊,通過概率統(tǒng)計度量給定長度窗口內(nèi)的模體組合的統(tǒng)計顯著性,相應(yīng)方法如MSCAN[22]和MCAST[28]等,或者使用組合方法,在指定窗口大小范圍內(nèi)搜索在多個序列共同出現(xiàn)一組模體實例的最小區(qū)域,將其作為候選順式調(diào)控模塊,如CMStalker[29]等.這類方法本質(zhì)上假定了序列窗內(nèi)的模體之間獨立同分布.這類方法雖然簡單直接,但需要合理確定窗大小以及度量統(tǒng)計顯著性的打分閾值,而這些參數(shù)在實際應(yīng)用中通常很難確定;此外,這類方法也忽略了順式調(diào)控模塊可能的調(diào)控結(jié)構(gòu)(或調(diào)控語法),如模體間的順序和距離.

另一類方法基于概率模型,通過對序列或順式調(diào)控模塊建立概率模型,進而找出待搜索的目標(biāo)序列中的順式調(diào)控模塊.基于概率模型的方法,除了少數(shù)采用判別模型的方法,如HexDiff[30]、Regulatory Potential[31]等外,大部分的方法使用生成模型,主要是隱馬爾科夫模型(Hidder Markov model, HMM).HMM模型的主要優(yōu)勢在于,它可以對順式調(diào)控模塊的出現(xiàn)進行可靠的統(tǒng)計度量,并能刻畫順式調(diào)控模塊的調(diào)控語法.此外,HMM模型所使用的期望最大化的參數(shù)估計算法,可以自動調(diào)節(jié)大量的參數(shù),避免了手動設(shè)置的麻煩.基于HMM的順式調(diào)控模塊識別方法通常將順式調(diào)控模塊看作由一組過表達的模體和背景組合生成的序列片段.與窗聚集方法相比,它們不僅考慮了構(gòu)成順式調(diào)控模塊的模體組合,也同時考慮了構(gòu)成順式調(diào)控模塊的模體之間的距離.最初的一些方法,如CisModule[20],使用HMM間接捕捉順式調(diào)控模塊內(nèi)部以及順式調(diào)控模塊之間背景的概率分布.但這種方法僅使用了一般的順式調(diào)控模塊內(nèi)部背景,并未推斷順式調(diào)控模塊內(nèi)的模體之間的任何順序.后續(xù)的方法,進一步擴展這種模型的表示,如Stubb[32]方法,創(chuàng)建了一個僅包含模體和背景兩個狀態(tài)的HMM模型,使用統(tǒng)計方法度量模體對間共同出現(xiàn)的顯著性,進而決定是否引入相應(yīng)的轉(zhuǎn)移概率.通過使用定義在指定長度窗口內(nèi)序列上的打分函數(shù)度量窗口內(nèi)模體聚集顯著性來預(yù)測順式調(diào)控模塊.該方法僅利用了HMM模型的轉(zhuǎn)移概率,沒有使用任何其他HMM模型特性.后來的方法CORECLUST[33]和TSHAS[14]建立了更復(fù)雜的HMM模型,引入了順式調(diào)控模塊內(nèi)部背景狀態(tài),加入了模體到模體的概率轉(zhuǎn)移,更細(xì)致地刻畫了順式調(diào)控模塊的調(diào)控結(jié)構(gòu).訓(xùn)練和解碼算法使用標(biāo)準(zhǔn)的Baum-Welch算法[34]和Veterbi算法[34].這類方法使用了HMM的特性,但局限于HMM的表達能力,建立的調(diào)控模型并不直觀,添加了大量的輔助狀態(tài).另一類HMM相關(guān)的順式調(diào)控模塊識別方法,使用加強的HMM來刻畫順式調(diào)控模塊的調(diào)控結(jié)構(gòu);如BayCis算法[35],使用貝葉斯層次HMM模型,對順式調(diào)控模塊和包含順式調(diào)控模塊調(diào)控序列進行建模.利用了層次HMM的特性建立了模體之間的概率轉(zhuǎn)移.該模型將HMM狀態(tài)轉(zhuǎn)移參數(shù)看作隨機變量,并引入貝葉斯先驗.該模型雖然結(jié)構(gòu)直觀,表達能力強,但模型的訓(xùn)練和解碼需要大量的計算.Lemnian等提出的基于Extended sunflower HMM的順式調(diào)控模塊識別方法[19],使用Extended sunflower HMM,對模型的刻畫深入到模體內(nèi)部,不僅刻畫了模體間的依賴,還刻畫了模體內(nèi)部堿基之間的依賴關(guān)系,但該方法僅用于同型順式調(diào)控模塊的識別.

此外,還有一類方法利用相近物種的進化保守性來識別順式調(diào)控模塊,如 MorphMS[26]、MultiModule[36]和ReLA[27]等.這類方法首先通過雙序列或多序列比對同源基因的調(diào)控區(qū)找出其中的保守區(qū)域,然后使用其他方法在保守區(qū)域中搜索順式調(diào)控模塊.由于大多數(shù)基因的調(diào)控區(qū)中存在著大量重復(fù)(Duplication)和改組(Shuffling)的序列片段,很難進行序列比對,所以這類方法并不總是有效.

將順式調(diào)控模塊的多種特征信息結(jié)合在一起,有助于提高識別順式調(diào)控模塊的準(zhǔn)確率.順式調(diào)控模塊雖然結(jié)構(gòu)具有不確定性,很難刻畫,但作為頻繁出現(xiàn)在多個被調(diào)控基因調(diào)控區(qū)中的調(diào)控功能單位,很可能包含某些保守成分.MOPAT[37]從一組同源基因的調(diào)控區(qū)中搜索保守模體模塊(這里的保守模體模塊被定義為在一組同源基因調(diào)控區(qū)中頻繁出現(xiàn)且具有一定距離約束的相鄰模體對,不考慮兩個模體的相對順序)出發(fā),查找包含多個這種保守模體模塊的區(qū)域,將其作為候選順式調(diào)控模塊,找出了一些保守的順式調(diào)控模塊.這一事實間接說明了順式調(diào)控模塊中存在保守成分.這里,我們同樣利用保守模體模塊這種順式調(diào)控模塊的保守成分,但限定保守模體模塊內(nèi)的模體具有特定的次序.然后,進一步將順式調(diào)控模塊保守性假設(shè)從同源基因(不同物種)推廣到共調(diào)控基因(同一物種).最終,我們將順式調(diào)控模塊表示為由單模體和保守模體模塊混合組成的,具有部分保守特征的調(diào)控結(jié)構(gòu)(也稱調(diào)控語法),從而將順式調(diào)控模塊結(jié)構(gòu)保守性和其內(nèi)部模體傾向于聚集的特征結(jié)合起來.為了刻畫這種復(fù)雜的順式調(diào)控模塊調(diào)控結(jié)構(gòu),我們使用一種被稱為Segmental HMM[38]的增強HMM模型來表達.

基于此,本文提出了一種識別順式調(diào)控模塊的概率模型方法SegHMC(Segmental HMM model for discovery of cis-regulatory module).該方法使用Segmental HMM在給定候選模體集上構(gòu)建同源或共調(diào)控基因調(diào)控區(qū)序列和順式調(diào)控模塊的調(diào)控語法結(jié)構(gòu).同一般的識別順式調(diào)控模塊的HMM模型相比,我們不僅將順式調(diào)控模塊表示為模體的組合,還將模體共同出現(xiàn)的頻率、模體順序偏好以及順式調(diào)控模塊中的相鄰模體之間距離分布等特征引入到順式調(diào)控模塊的調(diào)控語法當(dāng)中,這些特征可以有效提高順式調(diào)控模塊的識別精度.此外,為了處理真核生物基因長的調(diào)控區(qū),我們對模型進行了降低搜索空間的優(yōu)化.這種優(yōu)化通過提前進行片段分割,顯式建立Segmental HMM狀態(tài)轉(zhuǎn)換圖,去除了大量不必要的搜索路徑,降低了搜索空間,同時又不失精度.得到的模型可用于待搜索目標(biāo)基因調(diào)控區(qū)中甚至整個基因組中的相似順式調(diào)控模塊識別.我們分別在一個模擬數(shù)據(jù)集和兩個真實生物數(shù)據(jù)集: Muscle數(shù)據(jù)集和果蠅早期發(fā)育數(shù)據(jù)集上對我們的方法進行測試,并選取當(dāng)前主要方法進行比較,所有方法識別順式調(diào)控模塊的準(zhǔn)確率使用通用評價指標(biāo)相關(guān)系數(shù)(Correlation coefficient,CC)和F1-score來度量.實驗結(jié)果表明,我們的方法識別順式調(diào)控模塊的準(zhǔn)確率顯著優(yōu)于當(dāng)前的主要方法.

1 SegHMC算法

1.1 SegHMC的Segmental HMM模型

Segmental HMM[38]是HMM的一個擴展,也稱Generalized HMM.與一般HMM的每個狀態(tài)僅能發(fā)射一個堿基相比,Segmental HMM的每個狀態(tài)可以發(fā)射可變長度的堿基序列片段;狀態(tài)所發(fā)射的堿基序列可由一個片段模型來表示.該片段模型,給出了生成長度為u的觀察序列o=o1o2···ou的聯(lián)合概率,可由下式表示.

因此,片段模型由兩個分布組成,一個是描述片段長度似然的片段長度分布ds(u),另一個為表示不同長度觀察序列發(fā)射概率的發(fā)射模型es(o).因此,在順式調(diào)控模塊的識別模型中,可以根據(jù)對順式調(diào)控模塊和調(diào)控序列結(jié)構(gòu)的抽象,對這兩個分布給出具體的定義.

本文使用Segmental HMM,在片段層次上對順式調(diào)控模塊和調(diào)控序列的調(diào)控結(jié)構(gòu)進行建模,具有更強的表達能力;例如,可以對片段之間的依賴進行建模.本節(jié)將詳細(xì)闡述Segmental HMM模型,該模型主要包括:模型的構(gòu)建、狀態(tài)轉(zhuǎn)移概率、片段長度分布和生成狀態(tài)的發(fā)射模型.

1.1.1 Segmental HMM模型構(gòu)建

我們將轉(zhuǎn)錄調(diào)控序列的調(diào)控結(jié)構(gòu)定義如下.轉(zhuǎn)錄調(diào)控序列由一系列的順式調(diào)控模塊和順式調(diào)控模塊之間的背景(稱為全局背景)構(gòu)成,而每個順式調(diào)控模塊又由一組具有特定次序的模體和模體之間的背景(稱為局部背景)構(gòu)成,這種抽象具有明顯的層次性.基于這種結(jié)構(gòu)定義,給定的轉(zhuǎn)錄調(diào)控序列可由下列過程生成:

1)定位給定候選模體集中的模體在目標(biāo)轉(zhuǎn)錄調(diào)控序列中所有可能出現(xiàn)實例;

2)以這些被定位的模體實例為錨點,使用兩模體實例之間的背景(全局背景或局部背景,具體類別待定)序列連接這些模體實例,從而生成整個調(diào)控序列.

上述過程中,我們允許模體實例在空間上存在重疊,模體之間可能通過多種類型的背景序列相連;因此,存在許多平行的生成路徑.從這些路徑中,找出最可能的生成路徑,即可得出該轉(zhuǎn)錄調(diào)控序列最可能的調(diào)控結(jié)構(gòu),從而找出相應(yīng)的順式調(diào)控模塊.

將每個具體片段(模體、全局背景和局部背景)表示為Segmental HMM的一個狀態(tài),片段之間的連接對應(yīng)了兩個狀態(tài)之間的轉(zhuǎn)移,根據(jù)上述生成過程,我們顯式構(gòu)造Segmental HMM的狀態(tài)轉(zhuǎn)換圖.顯式構(gòu)建狀態(tài)轉(zhuǎn)換圖一方面移除了不必要的狀態(tài)路徑,減小算法的搜索空間;另一方面,更便于構(gòu)建模體的二元語法,模型順式調(diào)控模塊內(nèi)的相鄰模體間的一階依賴關(guān)系.為了標(biāo)識順式調(diào)控模塊,我們增加相應(yīng)的輔助狀態(tài):順式調(diào)控模塊開始狀態(tài)和順式調(diào)控模塊結(jié)束狀態(tài).圖2給出了表示一個調(diào)控序列調(diào)控語法結(jié)構(gòu)的Segmental HMM狀態(tài)轉(zhuǎn)換圖的具體例子,整個模型所包含的狀態(tài)如下:

2)模體狀態(tài)M={m1,m2,···,mK};

4)順式調(diào)控模塊狀態(tài)C,又由順式調(diào)控模塊開始狀態(tài)Cs和順式調(diào)控模塊結(jié)束狀態(tài)Ce構(gòu)成,即C=Cs∪Ce=

Segmental HMM狀態(tài)轉(zhuǎn)換圖的具體構(gòu)造過程如算法1所示.

算法1.Segmental HMM狀態(tài)轉(zhuǎn)換圖的構(gòu)造

輸入:一組Motif的PWM集PWMS,用于搜索Motif的p-value閾值和一個調(diào)控序列

輸出:狀態(tài)轉(zhuǎn)換圖的狀態(tài)集Q和這些狀態(tài)之間的連接集T

1)創(chuàng)建該模型的初始狀態(tài) S和終止?fàn)顟B(tài) E

2)對PWMS中每個PWM,在所給調(diào)控序列中找出小于給定p-value閾值所有Motif匹配

3)根據(jù)找出的Motif匹配對所給調(diào)控序列進行分割,標(biāo)記狀態(tài)類型,創(chuàng)建相應(yīng)的狀態(tài)集M,Bg和Cs

5)以狀態(tài)在序列中的位置為關(guān)鍵字對狀態(tài)集Q進行排序

6)Ce←?

7)Bc←?

8)for每個狀態(tài)qi∈Q do

9)if qi為模型的初始狀態(tài)

10)從Q中順序取出下一狀態(tài)qi+1

11)T←T∪{qi→qi+1}

12)else if qi為一個全局背景狀態(tài)then

13)從Q中找出qi位置之后的第一個全局背景狀態(tài),記為qj

14)T←T∪{qi→qj}

15)從Q中找出qi位置之后的第一個CRM初始狀態(tài),記為qj

16)T←T∪{qj→qi}

17)for qi的每個前端模體狀態(tài)m do

18)創(chuàng)建一個CRM終止?fàn)顟B(tài)ce

19)Ce←Ce∪{ce}

20)T←T∪{m→ce}

21)T←T∪{ce→qi}

22)else if qi是一個CRM初始狀態(tài)then

23)T←T∪{qi→mi}

24)else if qi是一個Motif狀態(tài)then

25)從Q中找出qi位置之后且不與它重疊的下一Motif狀態(tài),記為qj

26)創(chuàng)建一個局部背景狀態(tài)bc

27)Bc←Bc∪{bc}

28)T←T∪{qi→bc}

29)T←T∪{bc→qj}

30)Q←Q∪Ce∪Bc

31)以狀態(tài)在序列中的位置為關(guān)鍵字對狀態(tài)集Q進行重新排序

32)return Q和T

圖2 Segmental HMM狀態(tài)轉(zhuǎn)移圖Fig.2 The state transition diagram of segmental HMM

在我們的模型中使用位置權(quán)重矩陣(Position weight matrice,PWM)[39]表示相應(yīng)的模體,在算法中提前給定待搜索順式調(diào)控模塊所包含的可能模體集,所對應(yīng)的PWM集表示為PWMS.算法所要求的其他輸入包括:搜索模體的p-value,以及待建模的轉(zhuǎn)錄調(diào)控序列.

在上述Segmental HMM狀態(tài)轉(zhuǎn)換圖的構(gòu)造算法中,第1行創(chuàng)建模型的初始狀態(tài)和結(jié)束狀態(tài);第2行,根據(jù)所給p-value找出給定模體集中模體及其反向互補的模體在序列中所有的出現(xiàn)實例.第3行,根據(jù)所找出的模體實例,分別構(gòu)建模體狀態(tài)、全局背景狀態(tài)和順式調(diào)控模塊開始狀態(tài),對應(yīng)的狀態(tài)集分別為M、Bg和Cs.第5行對所有狀態(tài)的集合進行排序.第6～7行,分別初始化順式調(diào)控模塊結(jié)束狀態(tài)集Ce和模體間局部背景狀態(tài)集Bc.第8～31行,確定各狀態(tài)間的轉(zhuǎn)移,具體為:對于模型的初始狀態(tài)(對應(yīng)于第9～11行),只需連接下一任意有效狀態(tài);對于全局背景狀態(tài),需要與下一全局背景狀態(tài)(對應(yīng)于第13～14行)、相鄰順式調(diào)控模塊初始狀態(tài)(對應(yīng)于第15～16行)和順式調(diào)控模塊終止?fàn)顟B(tài)(對應(yīng)于第17～21行)相連;對于順式調(diào)控模塊初始狀態(tài),則只需連接到對應(yīng)的模體狀態(tài),這對應(yīng)于第22～23行;對于模體狀態(tài),對應(yīng)于第24～29行,找出后續(xù)與當(dāng)前狀態(tài)不重疊的模體狀態(tài),然后創(chuàng)建相應(yīng)的局部背景狀態(tài),并依次連接這些狀態(tài).第30～31行,將新創(chuàng)建的順式調(diào)控模塊結(jié)束狀態(tài)和模體間局部背景狀態(tài)加入到總狀態(tài)集Q,并對Q重新排序.

關(guān)于上述構(gòu)造算法的幾點說明:

1)在第3行中,對每個模體,分別創(chuàng)建了位于模體前后的兩個可能的全局背景狀態(tài).所創(chuàng)建的全局背景狀態(tài),僅有一端位置是確定的(即模體的前端或模體的后端),在后面第13～14行的操作中,會進一步將這些半連接的全局背景片段連接起來,形成一個大的全局背景;

2)在第5行和第31行中,按照兩個關(guān)鍵字(位置、狀態(tài)的類型)將前面分別生成的、無次序的各種狀態(tài)按照生成轉(zhuǎn)錄調(diào)控序列的時空順序進行排序,以便于后面確定各狀態(tài)的連接轉(zhuǎn)移的操作;

3)在算法中,我們顯式地創(chuàng)建順式調(diào)控模塊的開始狀態(tài)和結(jié)束狀態(tài),一方面,可以使結(jié)構(gòu)更清晰,另一方面,也便于在推斷時確定順式調(diào)控模塊的邊界;

4)對于第27行,為了簡化順式調(diào)控模塊的結(jié)構(gòu)表示模型,我們假設(shè)順式調(diào)控模塊內(nèi)相鄰的模體間是非重疊的;

5)假定所給序列的長度為T,在整個算法中,耗時的操作主要集中在:模體的查找,最壞時間復(fù)雜度為O(KT),其中K為PWM的個數(shù);查找后繼狀態(tài)的操作,最壞時間復(fù)雜度為O(T2);對狀態(tài)集的排序,最壞時間復(fù)雜度為O(T2).因此,整個算法的時間復(fù)雜度為O(T2).

1.1.2 狀態(tài)轉(zhuǎn)移概率

在Segmental HMM狀態(tài)轉(zhuǎn)換圖中,每個狀態(tài)對應(yīng)模體、全局背景或局部背景這些類型的一個實例,狀態(tài)之間的轉(zhuǎn)移概率即為相應(yīng)狀態(tài)類型之間轉(zhuǎn)移的概率.

對于順式調(diào)控模塊內(nèi)的模體狀態(tài)mi和模體狀態(tài)mj之間的轉(zhuǎn)移概率,可由下式估計得到:

其中,t(m)表示模體狀態(tài)(對應(yīng)于模體實例)m所對應(yīng)的模體類型,A為模體狀態(tài)間的轉(zhuǎn)移計數(shù).

對于全局背景類型狀態(tài)bg,其反映了在長的序列區(qū)域中出現(xiàn)順式調(diào)控模塊的概率.由狀態(tài)bg到順式調(diào)控模塊狀態(tài)的轉(zhuǎn)移概率,或順式調(diào)控模塊到bg狀態(tài)的轉(zhuǎn)移概率,由于即使在很長的調(diào)控序列中順式調(diào)控模塊的數(shù)目也相對較少,所能獲得的數(shù)據(jù)難以訓(xùn)練出可靠的模型參數(shù).為避免過擬合,可由經(jīng)驗估計得出,作為常量參數(shù),在系統(tǒng)運行時設(shè)定.

1.1.3 片段長度分布

全局背景長度和局部背景長度分別表示了順式調(diào)控模塊之間以及順式調(diào)控模塊內(nèi)的模體之間的空白區(qū)域的長度分布.對于全局背景狀態(tài)bg和局部背景狀態(tài)bc,我們假定其序列長度分別滿足期望為wg和wc的幾何分布;這種假定一方面反映了我們對順式調(diào)控模塊結(jié)構(gòu)的不確定性,另一方面又為模型順式調(diào)控模塊內(nèi)的模體的二元語法特征提供足夠的適應(yīng)性.在該假設(shè)下,背景序列長度為d的概率為:

這里,b表示bg或bc,wi表示wg或wc.

對于模體狀態(tài)m,由于模體所對應(yīng)的位置權(quán)重矩陣[39]是直接從數(shù)據(jù)庫中獲取的,其長度wm及其特定位置堿基的概率都是已知的,所以模體狀態(tài)上的序列長度d′的概率分布是特定的,即:

1.1.4 生成狀態(tài)的發(fā)射模型

在本文模型中,只有模體狀態(tài)、全局背景狀態(tài)和局部背景狀態(tài)為生成狀態(tài).每種生成狀態(tài)發(fā)射長度滿足特定分布的堿基序列片段.

對于全局背景和局部背景狀態(tài),我們分別使用k階Markov模型和m階的局部Markov模型.在局部Markov模型中,位置t處堿基的條件概率僅由以位置t為中心長度為2D的窗口內(nèi)的序列片段來估計.這可采用記筆記的方式預(yù)先計算出每個位置堿基的條件概率,并存儲計算的結(jié)果,需要時直接查表即得.

對于模體的生成概率,在本文模型中,使用經(jīng)典的PM 模型[40].假定模體實例為O,該模體的PWM Θ=[θ1,θ2,···,θL],其中θi(1≤i≤L)為堿基頻率的列向量,則模體狀態(tài)所對應(yīng)的堿基序列片段為O的概率為:

這里oi為模體實例O中第i位置的堿基.

1.2 解碼和訓(xùn)練算法

在我們的模型中,將輸入的序列分為訓(xùn)練集和測試集.在訓(xùn)練集上訓(xùn)練出模型參數(shù)后,使用已訓(xùn)練的模型識別給定測試集中所有序列的順式調(diào)控模塊,這一過程表現(xiàn)為解碼出模型的最優(yōu)狀態(tài)路徑過程.在Segmental HMM模型中,最優(yōu)狀態(tài)路徑可形式化定義為:

給定長度為T的轉(zhuǎn)錄調(diào)控序列(即觀測序列) O=o1o2···oT,設(shè)其對應(yīng)的狀態(tài)序列為Π=(π1,···,πT),則該轉(zhuǎn)錄調(diào)控序列所對應(yīng)的最優(yōu)狀態(tài)序列可表示為:

這里 Psi(di)表示狀態(tài) si的序列片段長度的概率分布,asi,si+1為狀態(tài)si到狀態(tài)si+1的轉(zhuǎn)移概率,e(oti+1···ti+1|si)表示狀態(tài)si生成觀測序列片段oti+1···ti+1的概率.

由于缺少足夠的標(biāo)注數(shù)據(jù),本文模型使用無監(jiān)督的Baum-Welch算法[34]直接從訓(xùn)練集中訓(xùn)練系統(tǒng)的模型參數(shù).對于模型的初始狀態(tài)概率,由于它只確定了輸入序列的第一個位置的初始功能狀態(tài),在沿序列的后續(xù)操作中,其影響完全可以忽略,在本文模型中簡單地由均勻分布隨機生成.

我們已提前標(biāo)出對應(yīng)片段的可能狀態(tài),創(chuàng)建了Segmental HMM的狀態(tài)轉(zhuǎn)換圖.因此,在解碼時,不再需要通過使用像最大似然之類的方法去推斷最可能的片段分割位置,可直接使用解碼算法找出最優(yōu)路徑.為了求式(6)所對應(yīng)的最優(yōu)狀態(tài)路徑,本文使用基于動態(tài)規(guī)劃的Veterbi算法[34],記為SegHMC Veterbi,并把它作為模型的缺省設(shè)置.此外,為了提供足夠的彈性,代替求最優(yōu)狀態(tài)路徑,本文還給出了類似于MAP(Maximum a posteriori probability)算法[34]基于閾值的后驗解碼算法,該算法給出了最可能的狀態(tài)路徑,記為SegHMC threshold.與MAP算法輸出每個后驗概率最大的序列區(qū)域相比,SegHMC threshold輸出后驗概率大于指定閾值包含順式調(diào)控模塊的序列區(qū)域.在SegHMC threshold算法中,本文搜索后驗概率大于給定閾值且至少包含兩個模體的連續(xù)區(qū)域作為候選順式調(diào)控模塊順式調(diào)控模塊.候選順式調(diào)控模塊區(qū)域的邊界定義為首個模體的起始位置,和最后一個模體的結(jié)束位置.在本文模型中,選擇的閾值范圍為[0.45,0.70],在模型的后驗推斷中該范圍內(nèi)的閾值通常能給出好的性能.相對于完全輸出后驗概率最大的MAP輸出來講,能通過合理地選取相應(yīng)的閾值在精度和召回率之間達到一個平衡.

1.3 SegHMC算法的整體框架

為了找出目標(biāo)序列中的順式調(diào)控模塊,本文算法需要給定相應(yīng)的輸入,它主要包括:訓(xùn)練集、測試集、常量參數(shù)、候選模體的PWM集、篩選順式調(diào)控模塊的閾值.算法通過執(zhí)行如下的過程,輸出所有測試集中大于給定閾值的順式調(diào)控模塊:

1)在訓(xùn)練集上,使用第1.1.1節(jié)中的Segmental HMM構(gòu)造算法構(gòu)造相應(yīng)的Segmental HMM模型,使用Baum-Welch算法訓(xùn)練出Segmental HMM的狀態(tài)轉(zhuǎn)移概率;

2)在測試集上,使用第1.1.1節(jié)中的Segmental HMM構(gòu)造算法構(gòu)造相應(yīng)的Segmental HMM模型,利用第1)步訓(xùn)練得出的模型參數(shù),使用Veterbi算法解碼或基于閾值的后驗解碼算法找出最優(yōu)或最可能的狀態(tài)路徑,預(yù)測出測試集中所有序列的順式調(diào)控模塊,進一步過濾得出大于給定閾值的順式調(diào)控模塊集.

模型的常量參數(shù)主要包括:搜索候選位點的pvalue;開始一個順式調(diào)控模塊的概率參數(shù)ps和結(jié)束一個順式調(diào)控模塊的概率參數(shù)pe;相鄰位點距離的幾何分布參數(shù)ml和相鄰順式調(diào)控模塊距離的幾何參數(shù)mg;以及順式調(diào)控模塊的權(quán)重閾值w.這些參數(shù)的取值首先結(jié)合真實生物數(shù)據(jù)特征的先驗知識來選取;例如,一個序列中通常的順式調(diào)控模塊的含量、順式調(diào)控模塊的平均長度、順式調(diào)控模塊內(nèi)模體間的平均距離等.其次,通過在模擬數(shù)據(jù)和搜集到的順式調(diào)控模塊數(shù)據(jù)集上,試驗一個取值范圍內(nèi)不同的取值,選擇在多數(shù)情況下能產(chǎn)生好的結(jié)果的值作為模型的參數(shù).

算法所輸出的順式調(diào)控模塊的具體信息主要包括:序列中所有找到的順式調(diào)控模塊的位置及其分值,以及構(gòu)成順式調(diào)控模塊的模體、相應(yīng)模體的位置和分值;其中,順式調(diào)控模塊分值和順式調(diào)控模塊內(nèi)的各個模體的分值,分別標(biāo)識了所找到的順式調(diào)控模塊和相應(yīng)模體的統(tǒng)計顯著性.我們將順式調(diào)控模塊的分值定義為對應(yīng)的序列片段由順式調(diào)控模塊狀態(tài)生成的后驗概率和該序列片段由全局背景狀態(tài)生成的后驗概率的log似然比值;類似地,模體的分值定義為該模體所對應(yīng)的序列片段由模體狀態(tài)生成的后驗概率和該序列片段由全局背景狀態(tài)生成的后驗概率的log似然比值.

2 實驗結(jié)果分析

2.1 實驗數(shù)據(jù)

我們在一個模擬數(shù)據(jù)集和兩個真實生物數(shù)據(jù)集上測試本文方法:脊椎動物的Muscle特定表達系統(tǒng)和果蠅早期胚胎發(fā)育系統(tǒng),以下簡稱Muscle數(shù)據(jù)集和果蠅早期發(fā)育數(shù)據(jù)集.這兩個真實生物數(shù)據(jù)集是評價當(dāng)前順式調(diào)控模塊方法性能的標(biāo)準(zhǔn)數(shù)據(jù)集[9?10],分別代表了兩類不同的序列數(shù)據(jù);其中Muscle數(shù)據(jù)上的序列為共調(diào)控序列,果蠅早期發(fā)育數(shù)據(jù)集中的序列為同源序列.

模擬數(shù)據(jù)集由30條序列和一個包含52個模體的模體集構(gòu)成,這52個模體從TRANSFAC數(shù)據(jù)庫中隨機提取;其中,每個序列的長度均為20kbp.序列集中的每條序列包含0～3個順式調(diào)控模塊;這些順式調(diào)控模塊的長度隨機為200～1500bp,模體位點間的空白的平均長度為50bp;每個順式調(diào)控模塊包含2～6種不同的模體,大約有15個模體位點實例;此外,為了模擬某些模體傾向于共同出現(xiàn)的偏好,設(shè)定每個順式調(diào)控模塊內(nèi)包含40%的有序模體對;順式調(diào)控模塊內(nèi)部和順式調(diào)控模塊間的背景序列均采用3階Markov模型,模型參數(shù)通過掃描D.melanogaster基因組間的序列獲得.

Muscle數(shù)據(jù)集最初由Wassermen和Fickett收集編錄[41];后來,Klepper等選取該數(shù)據(jù)集的一個子集并進行了擴展,作為文獻[9]中的基準(zhǔn)數(shù)據(jù).我們所使用的Muscle數(shù)據(jù)集即為文獻[9]中的基準(zhǔn)數(shù)據(jù),該數(shù)據(jù)集共包含5個模體和來自不同物種24條共調(diào)控序列.所包含的模體具體為Mef2、Myf、Sp1、SRF和Tef,這些模體在肌肉的轉(zhuǎn)錄調(diào)控中起重要作用.所包含的24條序列,平均長度為850bp,分別來自老鼠、人類、雞等物種.24條序列共包含5個模體的84個實例;每條序列均包含1個順式調(diào)控模塊,這些順式調(diào)控模塊的平均長度為120bp,其中最短的為14bp,最長的為294bp.

本文模型可以在一組具有相似調(diào)控結(jié)構(gòu)基因的調(diào)控區(qū)上進行訓(xùn)練,然后用于搜索其他基因甚至整個基因組中的相似順式調(diào)控模塊.但限于收集完整標(biāo)注全基因組數(shù)據(jù)的困難,我們僅在已標(biāo)注的調(diào)控果蠅早期前后軸發(fā)育系統(tǒng)的一組基因上進行了測試.相較于Muscle數(shù)據(jù)集,果蠅早期發(fā)育數(shù)據(jù)集中的基因有相對較長的順式調(diào)控模塊.在該數(shù)據(jù)集中,我們選取參與果蠅早期胚胎發(fā)育轉(zhuǎn)錄調(diào)控的7個模體:Bcd、Hb、Cad、Kr、Kni、Tll和Gt,并從iDMMPWM數(shù)據(jù)庫[42]中下載這些模體的PWM.我們從在果蠅早期胚胎發(fā)育過程中起重要作用的基因中選取一個包含7個基因的子集.這7個基因具體為:kni、kr、hb、tll、btd、eve和h,這些基因主要在果蠅早期胚胎的前后軸發(fā)育中調(diào)控果蠅幼體的體軸和體節(jié)發(fā)育.對于每個基因,我們選擇D.melanogaster和其12個同源的對應(yīng)基因序列作為搜索集.基因的染色體坐標(biāo)和同源信息從FlyBase數(shù)據(jù)庫[43]中獲取.典型的真核生物基因的調(diào)控區(qū)的長度通常為包含轉(zhuǎn)錄起始位點上游15kbps和下游5kbps的區(qū)域,但為了盡可能地覆蓋包含順式調(diào)控模塊的區(qū)域,我們將搜索區(qū)域定義為基因的轉(zhuǎn)錄起始位點附近的上游區(qū)20kbp和下游區(qū)20kbp,總共40kbp的序列區(qū)域,使用RepeatMasker軟件(http://www.repeatmasker.org)掩住(即替換為N字符串)搜索區(qū)內(nèi)的重復(fù)子序列.所選的7個模體和7個基因上的約91條序列構(gòu)成了我們實驗中所使用的果蠅早期發(fā)育數(shù)據(jù)集.我們從REDfly數(shù)據(jù)庫[44]中收集這些基因的順式調(diào)控模塊,并對存在重疊的順式調(diào)控模塊進行合并,將整理后的集合作為真實注釋的順式調(diào)控模塊基準(zhǔn)集.

2.2 評價

為了對本文方法在上述數(shù)據(jù)集上的預(yù)測結(jié)果進行客觀評價,并且使得評價結(jié)果不傾向于某一單個指標(biāo),我們使用評價指標(biāo)相關(guān)系數(shù)(Correlation coefficient,CC)[45]和定義在精度(Precision)和召回率(Recall)上的F1-score[46]作為堿基水平上的主要評價指標(biāo)對方法的預(yù)測準(zhǔn)確率進行評價,這些指標(biāo)被大多數(shù)順式調(diào)控識別方法所使用.此外,為了評價本文方法在位點水平上的精度/召回率,我們畫出了SegHMC threshold在不同閾值下的P/R曲線,以及模型缺省設(shè)置SegHMC Veterbi下的Precision/Recall值.這些評價指標(biāo)的具體定義為:

這里TP、FP、TN和FN分別表示真陽性、假陽性、真陰性和假陰性;式(9)中的精度(Precision)和召回率(Recall)的具體定義為:

1)精度(Precison):Precision=TP/(TP+ FP),度量了預(yù)測的正確率,表示預(yù)測結(jié)果中被正確識別的順式調(diào)控模塊的比率;

2)召回率(Recall):Recall=TP/(TP+FN),度量了順式調(diào)控模塊的識別率,表示基準(zhǔn)集中被正確識別的順式調(diào)控模塊的比率.

CC度量了預(yù)測集和基準(zhǔn)集中堿基位置的相關(guān)性,同時考慮了預(yù)測結(jié)果的真陽性率和真陰性率. CC取值范圍在–1和+1之間,+1表示預(yù)測結(jié)果與基準(zhǔn)集完全一致,–1表示預(yù)測結(jié)果與基準(zhǔn)集完全相反;當(dāng)預(yù)測結(jié)果接近隨機時,CC趨向于0.

從精度和召回率的定義,可以看到,它們分別衡量了預(yù)測算法的兩個方面:查全率和查準(zhǔn)率.雖然僅從定義上看這兩個指標(biāo)之間并沒有必然的聯(lián)系,但事實上它們在多數(shù)情況下是相互制約、相互矛盾的.這種制約和矛盾主要是由搜索策略的不完善造成的.通常情況下,如果希望獲得較高的召回率,則不得不降低搜索策略標(biāo)準(zhǔn),以盡可能地覆蓋基準(zhǔn)集中結(jié)果,從而引入了大量不相關(guān)的結(jié)果,降低了精度;如果想獲得較高的精度,則需要提高搜索策略標(biāo)準(zhǔn),這不可避免地過濾掉了相關(guān)的結(jié)果,造成了召回率的下降.因此,通常需要在這兩個指標(biāo)之間找到一個合理的平衡點.F1-score是精度和召回率的調(diào)和平均值,綜合了方法的精度和召回率,度量了方法在精度和召回率之間達到平衡的能力;其值在0和1之間,值越高,通常表示算法的性能越好(在真陽性率方面).

我們使用 Cython 語言,將本文方法SegHMC開發(fā)成相應(yīng)的工具,并在 Intel Xeon E5640@2.67GHz處理器,4GB內(nèi)存的Windows 7 64位系統(tǒng)的平臺上進行測試.

我們選取當(dāng)前的主要方法:BayCis[35]、Stubb[32]、MSCAN[22]、MotEvo[25]、ReLA[27]和CMStalker[29],在 Muscle 數(shù)據(jù)集上與本文方法SegHMC進行預(yù)測性能的比較.

在這6個被比較的方法中,BayCis和Stubb屬于概率模型方法,它們構(gòu)造相應(yīng)的HMM模型,訓(xùn)練模型參數(shù),使用訓(xùn)練的HMM解碼出給定序列中潛在的順式調(diào)控模塊.MSCAN是一個窗聚集方法,它計算給定序列滑動的窗內(nèi)出現(xiàn)的所有模體聚集的統(tǒng)計顯著性,將統(tǒng)計顯著性超過給定閾值的序列區(qū)域判定為候選順式調(diào)控模塊.MotEvo是窗聚集和概率模型的結(jié)合,它使用貝葉斯模型對滑動窗口內(nèi)的模體聚集區(qū)域進行打分,將分值大于給定閾值的區(qū)域作為候選順式調(diào)控模塊輸出.CMStalker是一個組合方法,它綜合了約束滿足規(guī)劃與參數(shù)松馳技術(shù),能有效遍歷模體組合的解空間,在給定的序列中找出可能的順式調(diào)控模塊.ReLA是一個基于局部比對的順式調(diào)控模塊搜索方法,它首先通過使用第三方模體搜索工具找出序列集中的所有模體,然后使用修改的Smith-Waterman算法,找出參考序列和被比較序列中得分最高的局部比對作為候選順式調(diào)控模塊.

2.3 模擬數(shù)據(jù)集上的結(jié)果

在模擬數(shù)據(jù)集上,模型固定參數(shù)具體設(shè)置為: ps=0.001,pe=0.1,p-value=0.01,和 w=100. mg和ml分別被設(shè)置為500bp和50bp,表示一個序列內(nèi)的順式調(diào)控模塊間的平均距離為500bp,順式調(diào)控模塊內(nèi)相鄰模體間的平均距離為50bp.對于其他被比較的方法,我們使用其默認(rèn)設(shè)置.

圖3(a)給出了所有方法的CC值和F1-score值.從圖中可以看出,我們的方法在這兩個評價指標(biāo)上明顯優(yōu)于其他方法.圖3(b)給出了所有方法的P/R值.從圖中可以看出,除MSCAN外的其他方法都有很高的召回率(Recall值),即所有方法的預(yù)測結(jié)果大都覆蓋了所植入的順式調(diào)控模塊,但預(yù)測的精度(Precision值)表現(xiàn)出很大的差異.但我們方法的兩個版本在保持高召回率的同時,精度超過絕大部分的其他方法,在精度和召回率之間達到一個很好的平衡.

圖3 SegHMC threshold和SegHMC Veterbi在模擬數(shù)據(jù)集上的預(yù)測性能Fig.3 The prediction performances of SegHMC threshold and SegHMC Veterbi on the synthetic dataset

2.4 Muscle數(shù)據(jù)集上的結(jié)果

考慮到該數(shù)據(jù)集上的序列所包含的順式調(diào)控模塊,大多具有模體高度聚集但總體長度較短的特征,我們設(shè)置模型參數(shù)ml為20bp,同時保持其他參數(shù)值不變.

圖4(a)給出了所有方法在Muscle數(shù)據(jù)集上的CC和F1-score值.從圖中可以看出,本文方法SegHMC在這兩個指標(biāo)上均優(yōu)于現(xiàn)存的方法(這里為模型缺省設(shè)置SegHMC Veterbi下的實驗結(jié)果).

圖4 SegHMC threshold和SegHMC Veterbi在Muscle集上的預(yù)測性能Fig.4 The prediction performances of SegHMC threshold and SegHMC Veterbi on the muscle dataset

為了考查SegHMC在該數(shù)據(jù)集上的具體表現(xiàn),我們根據(jù)不同的閾值給出SegHMC threshold在位點水平上預(yù)測性能的P/R曲線.對于模型缺省設(shè)置SegHMC Veterbi和其他被比較的方法,我們根據(jù)它們的缺省輸出在圖中給出相應(yīng)P/R值的點,如圖4(b)所示.在圖4(b)中,SegHMC threshold給出了一個平衡范圍,SegHMC Veterbi輸出位于曲線中間的一個點.

從圖4(b)中可以看出,本文方法的兩個版本在該數(shù)據(jù)集上較其他方法在預(yù)測精度和召回率之間均達到一個更好的平衡,具體為:在預(yù)測精度(Precision值)上高于其他大多數(shù)被比較的方法,召回率(Recall值)超過一半的方法.對于其他方法,分析圖中的數(shù)據(jù),還可以看到,在該數(shù)據(jù)集上不同的方法達到不同的P/R平衡,表現(xiàn)出不同的性能趨向.對于MotEvo和Stubb方法,它們傾向于給出較高召回率的預(yù)測,但整體上并未達到合理的平衡點,表現(xiàn)為預(yù)測精度(Precision值)明顯低于所有方法的平均值,僅約為本文方法的一半.很顯然,這種高的召回率是通過放寬順式調(diào)控模塊的篩選條件,盡可能覆蓋可能存在順式調(diào)控模塊的序列區(qū),以給出盡可能多的預(yù)測結(jié)果,犧牲預(yù)測精度換來的.而對于BayCis和CMStalker方法,為確保有較高的預(yù)測精度,它傾向于給出保守的預(yù)測結(jié)果,具體表現(xiàn)為其召回率(Recall值)明顯低于其他方法.對于MSCAN方法,它并未完全傾向于某一單個指標(biāo),而是在保證召回率的同時,給出了較高的預(yù)測精度,表現(xiàn)出良好的預(yù)測性能.在所有的方法中,ReLA并不傾向任何單個指標(biāo),給出了最保守的預(yù)測.

2.5 果蠅早期發(fā)育數(shù)據(jù)集上的結(jié)果

我們將這些方法進一步在果蠅早期發(fā)育數(shù)據(jù)集上進行測試,模型的所有常量參數(shù)設(shè)置與模擬數(shù)據(jù)集上相同.其中,對于Stubb,我們選擇其多物種版StubbMS(Stubb中的一個模塊,在以下的描述中我們?nèi)杂洖镾tubb).StubbMS利用相近物種間的保守性來提高方法的預(yù)測性能,具體為:通過雙序列比對,找出兩物種間的保守區(qū)域,對保守區(qū)內(nèi)的片段進行打分,將分值高于給定閾值的保守片段作為候選順式調(diào)控模塊.所有方法在該數(shù)據(jù)集上的整個測試方式如下:

1)對于 BayCis和 SegHMC,使用和D.melanogaster同源的所有其他果蠅物種的基因作為訓(xùn)練集,將D.melanogaster的對應(yīng)基因作為測試集(REDfly數(shù)據(jù)庫僅搜集了D.melanogaster基因的順式調(diào)控模塊標(biāo)注信息,其他同源物種的順式調(diào)控模塊暫時無法搜集);

2)對于需要雙序列比對的 Stubb方法,選取其原文中所使用的物種 (D.melanogaster和D.pseudoobscura),在相應(yīng)基因上進行測試;

3)對于僅在單序列上預(yù)測順式調(diào)控模塊的方法MSCAN和MotEvo,僅在D.melanogaster相應(yīng)的基因上進行測試;

4)對于ReLA,選取D.melanogaster相應(yīng)的基因作為參考序列,其余同源基因序列作為被比對的序列;

5)對于CMStalker,在所有同源基因上進行測試.

圖5(a)給出了所有方法在該數(shù)據(jù)集堿基水平上的CC和F1-score.從圖中可以看出,與其他方法相比,本文方法SegHMC在整個數(shù)據(jù)集的CC和F1-score值(這里為模型缺省設(shè)置SegHMC Veterbi下的實驗結(jié)果)上仍高于其他方法,表現(xiàn)出穩(wěn)定的預(yù)測性能.此外,通過對比可以發(fā)現(xiàn),與Muscle數(shù)據(jù)集上的結(jié)果相比,這些指標(biāo)的分值均有較大幅度的下降.考查這兩個數(shù)據(jù)集上數(shù)據(jù)的特征,發(fā)現(xiàn)這主要是由于搜索區(qū)長度的增加,而造成的信噪比(順式調(diào)控模塊相對于背景)的降低造成的.與Muscle數(shù)據(jù)集相比(序列平均長度為850bp),果蠅早期發(fā)育數(shù)據(jù)集中的序列更長(長度均為40kbp),各方法在識別單個模體時,不可避免地給出大量假陽性的預(yù)測,進而影響到順式調(diào)控模塊的預(yù)測結(jié)果,造成了預(yù)測結(jié)果中的大量假陽性預(yù)測.

圖5(b)給出了所有方法位點水平上的P/R曲線.從圖中可以看出,本文方法的兩個版本都達到了很好的平衡.并且閾值方法在某些閾值上較缺省設(shè)置達到更好的平衡,盡管如此,我們也應(yīng)該明白:在實際應(yīng)用中,我們所能達到的平衡完全基于數(shù)據(jù)的輸入,而P/R是完全未知的,我們很難選擇合適的閾值參數(shù),得到好的預(yù)測結(jié)果.因此,方法自動選擇合適的參數(shù)達到適當(dāng)?shù)腜/R平衡就顯得更為重要.從圖4(b)中,可以發(fā)現(xiàn):在Muslce數(shù)據(jù)集上表現(xiàn)良好的MSCAN方法,在該數(shù)據(jù)集上表現(xiàn)得很差,其精度(Precision值)僅處于中間水平;召回率(Recall值)明顯低于其他方法.在所有方法中MSCAN是唯一的基于窗聚集的方法,對于包含長度多變的順式調(diào)控模塊的序列,窗聚集方法很難確定一個合理的窗大小和打分閾值,從而給出好的預(yù)測結(jié)果(這里使用其默認(rèn)輸出);另一方面,MotEvo也是基于窗聚集的方法,但與純粹的基于窗聚集的MSCAN方法相比,它結(jié)合了貝葉斯概率模型對滑動窗內(nèi)的模體進行打分,從而有更好的性能,這也在某種程度上說明了,概率模型方法對不同類型的數(shù)據(jù)有更好的適應(yīng)性.對于引入系統(tǒng)發(fā)生的方法Stubb和ReLA,通過序列比對利用物種的保守性確實提高了預(yù)測精度(具有僅次于本文方法的Precision值),但完全基于物種間的保守性很難保證有較高的召回率(其Recall值最低).雖然我們的方法未直接通過序列比對利用物種間的保守性,但通過概率模型刻畫共調(diào)控或同源序列間結(jié)構(gòu)的相似性,在某種程度上也利用了物種間的系統(tǒng)發(fā)生關(guān)系,這也是本文方法優(yōu)于其他方法的一個原因.對于基于概率模型的方法BayCis,雖然在該數(shù)據(jù)集上給出了較高召回率(Recall值)的預(yù)測,但其預(yù)測的精度是最低的(Precision值).對于純粹組合搜索方法CMStalker,給出了很保守的預(yù)測,與在Muscle數(shù)據(jù)集上的結(jié)果相比,預(yù)測精度和召回率都下降很多.考慮到該數(shù)據(jù)集上的順式調(diào)控模塊的長度特點,我們推測可能超出了CMStalker的所能處理的空間范圍.

圖5 SegHMC threshold和SegHMC Veterbi在果蠅早期發(fā)育數(shù)據(jù)集上的預(yù)測性能Fig.5 The prediction performances of SegHMC threshold and SegHMC Veterbi on the early drosophila development dataset

為了考查不同的方法在該數(shù)據(jù)集中單個數(shù)據(jù)上的預(yù)測性能的差異,我們畫出了所有方法在該數(shù)據(jù)集中每個基因上的預(yù)測性能的箱線圖,如圖6所示.

圖6給出了每個方法CC值的中位數(shù)、下四分位數(shù)和上四分位數(shù)以及最小值和最大值.從圖中可以看出:所有方法在不同基因上的預(yù)測性能變化很大,這種變化要比方法總體性能之間的差異要大得多;這也說明了即使表現(xiàn)最好的方法在許多情況下也可能給出錯誤的預(yù)測.本文方法總體上較其他方法有更高的CC值,表現(xiàn)出更好的預(yù)測性能.雖然CMStalker在所有基因上的CC值變化最小,表現(xiàn)得最穩(wěn)定,但其CC值總體上明顯要低于大多數(shù)的方法.

圖6 所有方法在果蠅早期發(fā)育數(shù)據(jù)集中各基因上CC的變化Fig.6 A boxplot describing variation for all methods in CC across the genes in the early drosophila development dataset

2.6 引入的保守結(jié)構(gòu)信息帶來的性能提升

為了考查所引入的結(jié)構(gòu)信息(相鄰位點的相關(guān)性和位點間的距離分布)對方法性能的影響,我們回退模型到不考慮這些結(jié)構(gòu)信息一般模型,記為SegHMC-simple,在所有3個數(shù)據(jù)集上進行實驗對比,如圖7所示.

從圖中可以看出,加入結(jié)構(gòu)信息的模型SegHMC相對未加入相關(guān)信息的SegHMC-simple總體上具有明顯的性能提升,盡管對于不同的數(shù)據(jù)集性能提升不同.總體上,SegHMC在模擬數(shù)據(jù)和Muscle數(shù)據(jù)集上,較果蠅早期發(fā)育數(shù)據(jù)集上的提升要大.我們推測可能的原因在于,果蠅早期發(fā)育數(shù)據(jù)集上較大的搜索空間中所存在的大量噪聲造成的,即結(jié)構(gòu)信息的引入,可能增強原本的弱信號,引入了假陽性,抵消了一部分結(jié)構(gòu)信息所帶來的性能提升.檢查SegHMC在3個數(shù)據(jù)集上所找到的順式調(diào)控模塊,可以發(fā)現(xiàn),SegHMC能找出模擬數(shù)據(jù)集中所植入的70%以上共現(xiàn)模體對;在Muscle數(shù)據(jù)集中,找出了Mef2-Myf[7]、Myf-Sp1[7]以及Tef-Mef2[47]等經(jīng)過生物實驗驗證的模體對;在果蠅早期發(fā)育數(shù)據(jù)集上找出了共現(xiàn)的同型模體對[48],如Bcd-Bcd、Kr-Kr以及Hb-Hb等.

圖7 引入結(jié)構(gòu)信息后SegHMC在所有數(shù)據(jù)集上性能的提升Fig.7 An effect of inclusion of structural information on the prediction performance of SegHMC for all datasets

3 結(jié)論

本文將順式調(diào)控模塊的模體聚集特征和內(nèi)部結(jié)構(gòu)保守性結(jié)合起來,抽象出一種復(fù)雜的順式調(diào)控模塊的調(diào)控結(jié)構(gòu)表示(或稱調(diào)控語法),并借助于具有更強表達能力的Segmental HMM來表達,建立了一種識別順式調(diào)控模塊的概率模型.與其他基于概率類型的順式調(diào)控模塊識別方法相比,本文方法有如下的優(yōu)點:

1)我們不僅將順式調(diào)控模塊表示為模體的組合,還將模體共同出現(xiàn)的頻率、模體順序偏好以及順式調(diào)控模塊中的相鄰模體之間距離分布等特征引入到順式調(diào)控模塊的調(diào)控語法表示當(dāng)中,這些特征可以有效提高順式調(diào)控模塊的識別精度.

2)Segmental HMM的狀態(tài)可以表示一個序列片段,這一特性使得我們可以將所抽象的順式調(diào)控模塊的結(jié)構(gòu)元素(模體、背景等)直接映射為Segmental HMM的一個狀態(tài).在模體而不是堿基的抽象層次上對順式調(diào)控模塊的調(diào)控結(jié)構(gòu)進行建模,從而使得整個模型結(jié)構(gòu)的表示更清晰、更自然.此外,Segmental HMM并不限定片段長度的具體分布,使得我們可以根據(jù)模型假定選取特定的片段長度分布.

3)通過預(yù)先搜索模體實例,提前標(biāo)定相應(yīng)片段的類型,顯式構(gòu)建Segmental HMM的狀態(tài)轉(zhuǎn)換圖,有效減少了待搜索空間,提高了算法效率,使得本文方法可以處理真核生物普遍的長調(diào)控區(qū)序列.此外, segmental HMM模型所固有的在線(Online)性質(zhì),使得本文模型可以經(jīng)過在一組具有相似調(diào)控結(jié)構(gòu)基因的調(diào)控區(qū)上進行訓(xùn)練,然后用于搜索其他所要研究基因甚至整個基因組中的相似順式調(diào)控模塊.

為進一步提高識別順式調(diào)控模塊的性能,我們在后續(xù)的研究中打算從以下幾個方面著手,以進一步完善我們的工作:一方面,我們將搜集更多標(biāo)注的順式調(diào)控模塊,并使用更系統(tǒng)的方法,比如交叉驗證,來輔助模型參數(shù)的選擇.我們相信這些措施能使我們方法的性能得到進一步的提升.另一方面,當(dāng)前新的生物測序技術(shù),如染色質(zhì)共沉淀后微陣列分析(Chip-chip)或染色質(zhì)共沉淀后測序(Chip-seq)等的快速發(fā)展,產(chǎn)生了DNA雙螺旋結(jié)構(gòu)Profile、染色質(zhì)結(jié)構(gòu)、組蛋白修飾、蛋白質(zhì)占位等的大量實驗數(shù)據(jù),這些數(shù)據(jù)隱藏著基因表達調(diào)控規(guī)律,進一步分析提取這些可用于順式調(diào)控模塊預(yù)測的生物信息,并將這些信息與現(xiàn)有的計算模型結(jié)合,構(gòu)建更有效的順式調(diào)控模塊識別方法也是我們今后進一步努力的方向.

1 Wasserman W W,Sandelin A.Applied bioinformatics for the identification of regulatory elements.Nature Reviews Genetics,2004,5(4):276?287

2 Davidson E H.The Regulatory Genome:Gene Regulatory Networks in Development and Evolution.San Diego,California:Academic Press/Elsevier,2006.

3 Wang Pei,Lv Jin-Hu.Control of genetic regulatory networks: opportunities and challenges.Acta Automatica Sinica,2013,39(12):1969?1979 (王沛,呂金虎.基因調(diào)控網(wǎng)絡(luò)的控制:機遇與挑戰(zhàn).自動化學(xué)報, 2013,39(12):1969?1979)

4 Chen L N,Wang R S,Zhang X S.Biomolecular Networks: Methods and Applications in Systems Biology.Hoboken, New Jersey:Wiley,2009.

5 Kleinjan D A,Seawright A,Mella S,Carr C B,Tyas D A,Simpson T I,Mason J O,Price D J,van Heyningen V. Long-range downstream enhancers are essential for Pax6 expression.Developmental Biology,2006,299(2):563?581

6 Hardison R C,Taylor J.Genomic approaches towards finding cis-regulatory modules in animals.Nature Reviews Genetics,2012,13(7):469?483

7 Matys V,Kel-Margoulis O V,Fricke E,Liebich I,Land S,Barre-Dirrie A,Reuter I,Chekmenev D,Krull M,Hornischer K,Voss N,Stegmaier P,Lewicki-Potapov B,Saxel H,Kel A E,Wingender E.TRANSFAC?and its module TRANSCompel?:transcriptional gene regulation in eukaryotes.Nucleic Acids Research,2006,34(Database issue): D108?D110

8 Portales-Casamar E,Thongjuea S,Kwon A T,Arenillas D,Zhao X,Valen E,Yusuf D,Lenhard B,WassermanWW,Sandelin A.JASPAR 2010:the greatly expanded open-access database of transcription factor binding profiles.Nucleic Acids Research,2010,38(Database issue): D105?D110

9 Klepper K,Sandve G K,Abul O,Johansen J,Drablos F.Assessment of composite motif discovery methods.BMC Bioinformatics,2008,9:123

10 Su J,Teichmann S A,Down T A.Assessing computational methods of cis-regulatory module prediction.PLoS Computational Biology,2010,6(12):e1001020

12 Suryamohan K,Halfon M S.Identifying transcriptional cisregulatory modules in animal genomes.Wiley Interdisciplinary Reviews:Developmental Biology,2015,4(2):59?84

13 Thompson J A,Congdon C B.GAMI-CRM:using de novo motif inference to detect cis-regulatory modules.In:Proceedings of the 2014 IEEE Congress on Evolutionary Computation.Beijing,China:IEEE,2014.1022?1029

14 Zheng Shu-Rui.Research of Cis-regulatory Module Discovery Method Based on HMM Model[Master dissertation], Xidian University,China,2012 (鄭樹銳.基于HMM模型的順式調(diào)控模塊識別方法的研究[碩士學(xué)位論文],西安電子科技大學(xué),中國,2012)

15 Navarro C,Lopez F J,Cano C,Garcia-Alcalde F,Blanco A.CisMiner:genome-wide in-silico cis-regulatory module prediction by fuzzy itemset mining.PLoS One,2014,9(9): e108065

16 Rouault H,Santolini M,Schweisguth F,Hakim V.Imogene:identification of motifs and cis-regulatory modules underlying gene co-regulation.Nucleic Acids Research,2014, 42(10):6128?6145

17 Potier D,Seyres D,Guichard C,Iche-Torres M,Aerts S, Herrmann C,Perrin L.Identification of cis-regulatory modules encoding temporal dynamics during development.BMC Genomics,2014,15(1):534

18 Thompson J A,Congdon C B.Initial results in using de novo motif inference to detect cis-regulatory modules.In: Proceedings of the 2013 International Conference on Bioinformatics,Computational Biology and Biomedical Informatics.Washington DC,USA:ACM,2013.687

19 Lemnian I M,Eggeling R,Grosse I.Extended sunflower hidden Markov models for the recognition of homotypic cisregulatory modules.In:Proceedings of the 2013 German Conference on Bioinformatics.Gottingen,Germany,2013. 101?109

20 Zhou Q,Wong W H.CisModule:de novo discovery of cisregulatory modules by hierarchical mixture modeling.Proceedings of the National Academy of Sciences of the United States of America,2004,101(33):12114?12119

21 Gan Y L,Guan J H,Zhou S G,Zhang W X.Identifying cisregulatory elements and modules using conditional random fields.IEEE/ACM Transactions on Computational Biology and Bioinformatics,2014,11(1):73?82

22 Alkema W B,Johansson O,Lagergren J,Wasserman W W. MSCAN:identification of functional clusters of transcription factor binding sites.Nucleic Acids Research,2004,32(Web Server issue):W195?W198

23 Aerts S,Van Loo P,Thijs G,Moreau Y,De Moor B.Computational detection of cis-regulatory modules.Bioinformatics, 2003,19(Suppl 2):ii5?ii14

24 Sharan R,Ovcharenko I,Ben-Hur A,Karp R M.CREME:a framework for identifying cis-regulatory modules in humanmouse conserved segments.Bioinformatics,2003,19(Suppl 1):i283?i291

25 Arnold P,Erb I,Pachkov M,Molina N,van Nimwegen E. MotEvo:integrated Bayesian probabilistic methods for inferring regulatory sites and motifs on multiple alignments of DNA sequences.Bioinformatics,2012,28(4):487?494

26 Sinha S,He X.MORPH:probabilistic alignment combined with hidden Markov models of cis-regulatory modules.PLoS Computational Biology,2007,3(11):e216

28 Bailey T L,Noble W S.Searching for statistically significant regulatory modules.Bioinformatics,2003,19(Suppl 2):ii16?ii25

29 Leoncini M,Montangero M,Pellegrini M,Tillan K P.CMStalker:a combinatorial tool for composite motif discovery. IEEE/ACM Transactions on Computational Biology and Bioinformatics,2015,12(5):1123?1136

30 Chan B Y,Kibler D.Using hexamers to predict cisregulatory motifs in Drosophila.BMC Bioinformatics,2005, 6:262

31 Kolbe D,Taylor J,Elnitski L,Eswara P,Li J,Miller W, Hardison R,Chiaromonte F.Regulatory potential scores from genome-wide three-way alignments of human,mouse, and rat.Genome Research,2004,14(4):700?707

32 Sinha S,van Nimwegen E,Siggia E D.A probabilistic method to detect regulatory modules.Bioinformatics,2003, 19(Suppl 1):i292?i301

33 Nikulova A A,Favorov A V,Sutormin R A,Makeev V J, Mironov A A.CORECLUST:identification of the conserved CRM grammar together with prediction of gene regulation. Nucleic Acids Research,2012,40(12):e93

34 Durbin R,Eddy S R,Krogh A,Mitchison G.Biological Sequence Analysis:Probabilistic Models of Proteins and Nucleic Acids.Cambridge:Cambridge University Press,1998.

35 Lin T H,Ray P,Sandve G K,Uguroglu S,Xing E P.Bay-Cis:a Bayesian hierarchical HMM for cis-regulatory module decoding in metazoan genomes.In:Proceedings of the 12th Annual International Conference on Research in Computational Molecular Biology.Singapore:Springer,2008.66?81

36 Zhou Q,Wong W H.Coupling hidden Markov models for the discovery of Cis-regulatory modules in multiple species. Annals of Applied Statistics,2007,1(1):36?65

37 Hu J F,Hu H Y,Li X M.MOPAT:a graph-based method to predict recurrent cis-regulatory modules from known motifs. Nucleic Acids Research,2008,36(13):4488?4497

38 Russell M J.A segmental HMM for speech pattern modelling.In:Processing of the 1993 IEEE International Conference on Acoustics,Speech,and Signal Processing.Minneapolis,MN,USA:IEEE,1993.499?502

39 Stormo G D.DNA binding sites:representation and discovery.Bioinformatics,2000,16(1):16?23

40 Liu X,Brutlag D L,Liu J S.BioProspector:discovering conserved DNA motifs in upstream regulatory regions of co-expressed genes.In:Proceedings of the 6th Pacific Symposium on Biocomputing.Hawaii,USA,2001.127?138

41 Wasserman W W,Fickett J W.Identification of regulatory regions which confer muscle-specific gene expression.Journal of Molecular Biology,1998,278(1):167?181

42 Kulakovskiy I V,Makeev V J.Discovery of DNA motifs recognized by transcription factors through integration of different experimental sources.Biophysics,2009,54(6): 667?674

43 Tweedie S,Ashburner M,Falls K,Leyland P,McQuilton P, Marygold S,Millburn G,Osumi-Sutherland D,Schroeder A,Seal R,Zhang H,The FlyBase Consortium.FlyBase: enhancing Drosophila Gene Ontology annotations.Nucleic Acids Research,2009,37(Database issue):D555?D559

44 Gallo S M,Gerrard D T,Miner D,Simich M,Des Soye B, Bergman C M,Halfon M S.REDfly v3.0:toward a comprehensive database of transcriptional regulatory elements in Drosophila.Nucleic Acids Research,2011,39(Database issue):D118-D123

45 Tompa M,Li N,Bailey T L,Church G M,De Moor B,Eskin E,Favorov A V,Frith M C,Fu Y T,Kent W J,Makeev V J,Mironov A A,Noble W S,Pavesi G,Pesole G,Rgnier M,Simonis N,Sinha S,Thijs G,van Helden J,Vandenbogaert M,Weng Z P,Workman C,Ye C,Zhu Z.Assessing computational tools for the discovery of transcription factor binding sites.Nature Biotechnology,2005,23:137?144

46 Shaw W M Jr,Burgin R,Howell P.Performance standards and evaluations in IR test collections:cluster-based retrieval models.Information Processing&Management, 1997,33(1):1?14

47 Maeda T,Gupta M P,Stewart A F R.TEF-1 and MEF2 transcription factors interact to regulate muscle-specific promoters.Biochemical and Biophysical Research Communications,2002,294(4):791?797

48 Lifanov A P,Makeev V J,Nazina A G,Papatsenko D A. Homotypic regulatory clusters in Drosophila.Genome Research,2003,13(4):579?588

郭海濤 西安電子科技大學(xué)計算機學(xué)院博士研究生.主要研究方向為生物信息學(xué)算法,并行算法.

E-mail:ght_311@sina.com

(GUO Hai-Tao Ph.D.candidate at the School of Computer Science and Technology,Xidian University.His research interest covers bioinformatics algorithms and pararllel algorithms.)

霍紅衛(wèi) 博士,西安電子科技大學(xué)計算機學(xué)院教授.主要研究方向為大數(shù)據(jù)算法與壓縮數(shù)據(jù)結(jié)構(gòu),生物信息學(xué)算法,算法工程.本文通信作者.

E-mail:hwhuo@mail.xidian.edu.cn

(HUO Hong-Wei Ph.D.,professor at the School of Computer Science and Technology,Xidian University.Her research interest covers big data algorithms and compressed data structures,bioinformatics algorithms,algorithm engineering.Corresponding author of this paper.)

于 強 博士,西安電子科技大學(xué)計算機學(xué)院講師.主要研究方向為生物信息學(xué)算法,并行算法.

E-mail:qyu@mail.xidian.edu.cn

(YU Qiang Ph.D.,lecturer at the School of Computer Science and Technology,Xidian University. His research interest covers bioinformatics algorithms and pararllel algorithms.)

SegHMC:an Algorithm for Discovery of Cis-regulatory Module Based on Segmental HMM

GUO Hai-Tao1HUO Hong-Wei1YU Qiang1

Cis-regulatory module(CRM)plays a key role in metazoan gene transcriptional regulation,and the discovery of cis-regulatory module has been a crucial research topic recently.Many computational methods have been proposed to predict the cis-regulatory module,but it is still a main task to further improve the prediction accuracy for cis-regulatory modules.Combining multiple features of cis-regulatory module together can improve the prediction accuracy for cisregulatory module.Based on this,the paper presents an algorithm SegHMC(Segmental HMM model for discovery of cis-regulatory module)for the discovery of cis-regulatory module based on segmental HMM.The model further extends the representation of the structure of cis-regulatory module(or regulatory grammar),which not only describes a CRM as a combination of a group of motifs but also further introduces the frequency of the occurrence of motifs,the favour of the order of motifs,and the distance distribution between the adjacent motifs and other features.Experiments on the benchmark datasets demonstrate that the proposed algorithm outperforms the present main algorithms in the prediction accuracy.

Gene transcriptional regulation,motif,segmental HMM,discovery of cis-regulatory module

郭海濤,霍紅衛(wèi),于強.SegHMC:一種基于Segmental HMM模型的順式調(diào)控模塊識別算法.自動化學(xué)報,2016, 42(11):1718?1731

Guo Hai-Tao,Huo Hong-Wei,Yu Qiang.SegHMC:an algorithm for discovery of cis-regulatory module based on segmental HMM.Acta Automatica Sinica,2016,42(11):1718?1731

2015-05-18 錄用日期2016-06-06

Manuscript received May 18,2015;accepted June 6,2016

國家自然科學(xué)基金(61173025,61373044,61502366),中國博士后科學(xué)基金(2015M582621)資助

Supported by National Natural Science Foundation of China (61173025,61373044,61502366),the Chinese Postdoctoral Science Foundation(2015M582621)

本文責(zé)任編委黃慶明

Recommended by Associate Editor HUANG Qing-Ming

1.西安電子科技大學(xué)計算機學(xué)院西安710071

1.School of Computer Science and Technology,Xidian University,Xi'an 710071

DOI 10.16383/j.aas.2016.c150309