• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2016-12-17 08:23:58郭海濤霍紅衛(wèi)于強(qiáng)
    自動(dòng)化學(xué)報(bào) 2016年11期
    關(guān)鍵詞:模體集上調(diào)控

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

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

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

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

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

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

    圖1 順式調(diào)控模塊結(jié)構(gòu)示意圖(順式調(diào)控模塊是包含多個(gè)轉(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.)

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

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

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

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

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

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

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

    將順式調(diào)控模塊的多種特征信息結(jié)合在一起,有助于提高識(shí)別順式調(diào)控模塊的準(zhǔn)確率.順式調(diào)控模塊雖然結(jié)構(gòu)具有不確定性,很難刻畫,但作為頻繁出現(xiàn)在多個(gè)被調(diào)控基因調(diào)控區(qū)中的調(diào)控功能單位,很可能包含某些保守成分.MOPAT[37]從一組同源基因的調(diào)控區(qū)中搜索保守模體模塊(這里的保守模體模塊被定義為在一組同源基因調(diào)控區(qū)中頻繁出現(xiàn)且具有一定距離約束的相鄰模體對(duì),不考慮兩個(gè)模體的相對(duì)順序)出發(fā),查找包含多個(gè)這種保守模體模塊的區(qū)域,將其作為候選順式調(diào)控模塊,找出了一些保守的順式調(diào)控模塊.這一事實(shí)間接說明了順式調(diào)控模塊中存在保守成分.這里,我們同樣利用保守模體模塊這種順式調(diào)控模塊的保守成分,但限定保守模體模塊內(nèi)的模體具有特定的次序.然后,進(jìn)一步將順式調(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]的增強(qiáng)HMM模型來表達(dá).

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

    1 SegHMC算法

    1.1 SegHMC的Segmental HMM模型

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

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

    本文使用Segmental HMM,在片段層次上對(duì)順式調(diào)控模塊和調(diào)控序列的調(diào)控結(jié)構(gòu)進(jìn)行建模,具有更強(qiáng)的表達(dá)能力;例如,可以對(duì)片段之間的依賴進(jìn)行建模.本節(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)成,而每個(gè)順式調(diào)控模塊又由一組具有特定次序的模體和模體之間的背景(稱為局部背景)構(gòu)成,這種抽象具有明顯的層次性.基于這種結(jié)構(gòu)定義,給定的轉(zhuǎn)錄調(diào)控序列可由下列過程生成:

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

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

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

    將每個(gè)具體片段(模體、全局背景和局部背景)表示為Segmental HMM的一個(gè)狀態(tài),片段之間的連接對(duì)應(yīng)了兩個(gè)狀態(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)識(shí)順式調(diào)控模塊,我們?cè)黾酉鄳?yīng)的輔助狀態(tài):順式調(diào)控模塊開始狀態(tài)和順式調(diào)控模塊結(jié)束狀態(tài).圖2給出了表示一個(gè)調(diào)控序列調(diào)控語法結(jié)構(gòu)的Segmental HMM狀態(tài)轉(zhuǎn)換圖的具體例子,整個(gè)模型所包含的狀態(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閾值和一個(gè)調(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)對(duì)PWMS中每個(gè)PWM,在所給調(diào)控序列中找出小于給定p-value閾值所有Motif匹配

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

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

    6)Ce←?

    7)Bc←?

    8)for每個(gè)狀態(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為一個(gè)全局背景狀態(tài)then

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

    14)T←T∪{qi→qj}

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

    16)T←T∪{qj→qi}

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

    18)創(chuàng)建一個(gè)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是一個(gè)CRM初始狀態(tài)then

    23)T←T∪{qi→mi}

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

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

    26)創(chuàng)建一個(gè)局部背景狀態(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)鍵字對(duì)狀態(tài)集Q進(jìn)行重新排序

    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)控模塊所包含的可能模體集,所對(duì)應(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找出給定模體集中模體及其反向互補(bǔ)的模體在序列中所有的出現(xiàn)實(shí)例.第3行,根據(jù)所找出的模體實(shí)例,分別構(gòu)建模體狀態(tài)、全局背景狀態(tài)和順式調(diào)控模塊開始狀態(tài),對(duì)應(yīng)的狀態(tài)集分別為M、Bg和Cs.第5行對(duì)所有狀態(tài)的集合進(jìn)行排序.第6~7行,分別初始化順式調(diào)控模塊結(jié)束狀態(tài)集Ce和模體間局部背景狀態(tài)集Bc.第8~31行,確定各狀態(tài)間的轉(zhuǎn)移,具體為:對(duì)于模型的初始狀態(tài)(對(duì)應(yīng)于第9~11行),只需連接下一任意有效狀態(tài);對(duì)于全局背景狀態(tài),需要與下一全局背景狀態(tài)(對(duì)應(yīng)于第13~14行)、相鄰順式調(diào)控模塊初始狀態(tài)(對(duì)應(yīng)于第15~16行)和順式調(diào)控模塊終止?fàn)顟B(tài)(對(duì)應(yīng)于第17~21行)相連;對(duì)于順式調(diào)控模塊初始狀態(tài),則只需連接到對(duì)應(yīng)的模體狀態(tài),這對(duì)應(yīng)于第22~23行;對(duì)于模體狀態(tài),對(duì)應(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,并對(duì)Q重新排序.

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

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

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

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

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

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

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

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

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

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

    對(duì)于全局背景類型狀態(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ù)目也相對(duì)較少,所能獲得的數(shù)據(jù)難以訓(xùn)練出可靠的模型參數(shù).為避免過擬合,可由經(jīng)驗(yàn)估計(jì)得出,作為常量參數(shù),在系統(tǒng)運(yùn)行時(shí)設(shè)定.

    1.1.3 片段長度分布

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

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

    對(duì)于模體狀態(tài)m,由于模體所對(duì)應(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ā)射長度滿足特定分布的堿基序列片段.

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

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

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

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

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

    給定長度為T的轉(zhuǎn)錄調(diào)控序列(即觀測序列) O=o1o2···oT,設(shè)其對(duì)應(yīng)的狀態(tài)序列為Π=(π1,···,πT),則該轉(zhuǎn)錄調(diào)控序列所對(duì)應(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ù).對(duì)于模型的初始狀態(tài)概率,由于它只確定了輸入序列的第一個(gè)位置的初始功能狀態(tài),在沿序列的后續(xù)操作中,其影響完全可以忽略,在本文模型中簡單地由均勻分布隨機(jī)生成.

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

    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àn)解碼算法找出最優(yōu)或最可能的狀態(tài)路徑,預(yù)測出測試集中所有序列的順式調(diào)控模塊,進(jìn)一步過濾得出大于給定閾值的順式調(diào)控模塊集.

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

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

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

    2.1 實(shí)驗(yàn)數(shù)據(jù)

    我們?cè)谝粋€(gè)模擬數(shù)據(jù)集和兩個(gè)真實(shí)生物數(shù)據(jù)集上測試本文方法:脊椎動(dòng)物的Muscle特定表達(dá)系統(tǒng)和果蠅早期胚胎發(fā)育系統(tǒng),以下簡稱Muscle數(shù)據(jù)集和果蠅早期發(fā)育數(shù)據(jù)集.這兩個(gè)真實(shí)生物數(shù)據(jù)集是評(píng)價(jià)當(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條序列和一個(gè)包含52個(gè)模體的模體集構(gòu)成,這52個(gè)模體從TRANSFAC數(shù)據(jù)庫中隨機(jī)提取;其中,每個(gè)序列的長度均為20kbp.序列集中的每條序列包含0~3個(gè)順式調(diào)控模塊;這些順式調(diào)控模塊的長度隨機(jī)為200~1500bp,模體位點(diǎn)間的空白的平均長度為50bp;每個(gè)順式調(diào)控模塊包含2~6種不同的模體,大約有15個(gè)模體位點(diǎn)實(shí)例;此外,為了模擬某些模體傾向于共同出現(xiàn)的偏好,設(shè)定每個(gè)順式調(diào)控模塊內(nèi)包含40%的有序模體對(duì);順式調(diào)控模塊內(nèi)部和順式調(diào)控模塊間的背景序列均采用3階Markov模型,模型參數(shù)通過掃描D.melanogaster基因組間的序列獲得.

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

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

    2.2 評(píng)價(jià)

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

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

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

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

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

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

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

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

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

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

    圖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)控模塊,大多具有模體高度聚集但總體長度較短的特征,我們?cè)O(shè)置模型參數(shù)ml為20bp,同時(shí)保持其他參數(shù)值不變.

    圖4(a)給出了所有方法在Muscle數(shù)據(jù)集上的CC和F1-score值.從圖中可以看出,本文方法SegHMC在這兩個(gè)指標(biāo)上均優(yōu)于現(xiàn)存的方法(這里為模型缺省設(shè)置SegHMC Veterbi下的實(shí)驗(yàn)結(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在位點(diǎn)水平上預(yù)測性能的P/R曲線.對(duì)于模型缺省設(shè)置SegHMC Veterbi和其他被比較的方法,我們根據(jù)它們的缺省輸出在圖中給出相應(yīng)P/R值的點(diǎn),如圖4(b)所示.在圖4(b)中,SegHMC threshold給出了一個(gè)平衡范圍,SegHMC Veterbi輸出位于曲線中間的一個(gè)點(diǎn).

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

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

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

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

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

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

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

    5)對(duì)于CMStalker,在所有同源基因上進(jìn)行測試.

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

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

    圖6給出了每個(gè)方法CC值的中位數(shù)、下四分位數(shù)和上四分位數(shù)以及最小值和最大值.從圖中可以看出:所有方法在不同基因上的預(yù)測性能變化很大,這種變化要比方法總體性能之間的差異要大得多;這也說明了即使表現(xiàn)最好的方法在許多情況下也可能給出錯(cuò)誤的預(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)信息(相鄰位點(diǎn)的相關(guān)性和位點(diǎn)間的距離分布)對(duì)方法性能的影響,我們回退模型到不考慮這些結(jié)構(gòu)信息一般模型,記為SegHMC-simple,在所有3個(gè)數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn)對(duì)比,如圖7所示.

    從圖中可以看出,加入結(jié)構(gòu)信息的模型SegHMC相對(duì)未加入相關(guān)信息的SegHMC-simple總體上具有明顯的性能提升,盡管對(duì)于不同的數(shù)據(jù)集性能提升不同.總體上,SegHMC在模擬數(shù)據(jù)和Muscle數(shù)據(jù)集上,較果蠅早期發(fā)育數(shù)據(jù)集上的提升要大.我們推測可能的原因在于,果蠅早期發(fā)育數(shù)據(jù)集上較大的搜索空間中所存在的大量噪聲造成的,即結(jié)構(gòu)信息的引入,可能增強(qiáng)原本的弱信號(hào),引入了假陽性,抵消了一部分結(jié)構(gòu)信息所帶來的性能提升.檢查SegHMC在3個(gè)數(shù)據(jù)集上所找到的順式調(diào)控模塊,可以發(fā)現(xiàn),SegHMC能找出模擬數(shù)據(jù)集中所植入的70%以上共現(xiàn)模體對(duì);在Muscle數(shù)據(jù)集中,找出了Mef2-Myf[7]、Myf-Sp1[7]以及Tef-Mef2[47]等經(jīng)過生物實(shí)驗(yàn)驗(yàn)證的模體對(duì);在果蠅早期發(fā)育數(shù)據(jù)集上找出了共現(xiàn)的同型模體對(duì)[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)控語法),并借助于具有更強(qiáng)表達(dá)能力的Segmental HMM來表達(dá),建立了一種識(shí)別順式調(diào)控模塊的概率模型.與其他基于概率類型的順式調(diào)控模塊識(shí)別方法相比,本文方法有如下的優(yōu)點(diǎn):

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

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

    3)通過預(yù)先搜索模體實(shí)例,提前標(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ū)上進(jìn)行訓(xùn)練,然后用于搜索其他所要研究基因甚至整個(gè)基因組中的相似順式調(diào)控模塊.

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

    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ò)的控制:機(jī)遇與挑戰(zhàn).自動(dòng)化學(xué)報(bào), 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)控模塊識(shí)別方法的研究[碩士學(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é)計(jì)算機(jī)學(xué)院博士研究生.主要研究方向?yàn)樯镄畔W(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é)計(jì)算機(jī)學(xué)院教授.主要研究方向?yàn)榇髷?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.)

    于 強(qiáng) 博士,西安電子科技大學(xué)計(jì)算機(jī)學(xué)院講師.主要研究方向?yàn)樯镄畔W(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),于強(qiáng).SegHMC:一種基于Segmental HMM模型的順式調(diào)控模塊識(shí)別算法.自動(dòng)化學(xué)報(bào),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é)計(jì)算機(jī)學(xué)院西安710071

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

    DOI 10.16383/j.aas.2016.c150309

    猜你喜歡
    模體集上調(diào)控
    基于Matrix Profile的時(shí)間序列變長模體挖掘
    Cookie-Cutter集上的Gibbs測度
    鏈完備偏序集上廣義向量均衡問題解映射的保序性
    如何調(diào)控困意
    經(jīng)濟(jì)穩(wěn)中有進(jìn) 調(diào)控托而不舉
    中國外匯(2019年15期)2019-10-14 01:00:34
    植入(l, d)模體發(fā)現(xiàn)若干算法的實(shí)現(xiàn)與比較
    復(fù)扇形指標(biāo)集上的分布混沌
    基于網(wǎng)絡(luò)模體特征攻擊的網(wǎng)絡(luò)抗毀性研究
    順勢而導(dǎo) 靈活調(diào)控
    基于模體演化的時(shí)序鏈路預(yù)測方法
    亚洲av片天天在线观看| 自线自在国产av| 人妻 亚洲 视频| 天天添夜夜摸| 极品人妻少妇av视频| 两个人看的免费小视频| 亚洲熟女毛片儿| 国产不卡一卡二| 国产精品九九99| 十八禁人妻一区二区| 日韩中文字幕视频在线看片| 一区二区三区乱码不卡18| 中文字幕av电影在线播放| 正在播放国产对白刺激| 91成年电影在线观看| 人妻 亚洲 视频| 国产日韩欧美在线精品| 韩国精品一区二区三区| 侵犯人妻中文字幕一二三四区| 久久久久网色| 欧美激情久久久久久爽电影 | 香蕉久久夜色| 99久久99久久久精品蜜桃| 色播在线永久视频| 天天躁日日躁夜夜躁夜夜| 午夜免费鲁丝| 蜜桃国产av成人99| 亚洲熟妇熟女久久| av福利片在线| 伦理电影免费视频| 999久久久精品免费观看国产| 久久精品国产亚洲av香蕉五月 | 在线av久久热| 99国产精品一区二区三区| 999久久久国产精品视频| 深夜精品福利| 免费观看a级毛片全部| 黄片播放在线免费| 国产成人欧美| 亚洲av日韩精品久久久久久密| 成年人黄色毛片网站| 亚洲欧洲日产国产| 亚洲五月色婷婷综合| 精品视频人人做人人爽| 中国美女看黄片| svipshipincom国产片| 在线十欧美十亚洲十日本专区| 免费看a级黄色片| 视频区欧美日本亚洲| 欧美乱妇无乱码| 国产精品一区二区在线不卡| 妹子高潮喷水视频| 在线天堂中文资源库| 欧美精品一区二区免费开放| 国产免费视频播放在线视频| 十八禁高潮呻吟视频| 欧美日韩中文字幕国产精品一区二区三区 | 美女主播在线视频| 真人做人爱边吃奶动态| 亚洲美女黄片视频| 性色av乱码一区二区三区2| 性少妇av在线| 亚洲va日本ⅴa欧美va伊人久久| 国产男靠女视频免费网站| 国产免费现黄频在线看| 最新在线观看一区二区三区| 色婷婷久久久亚洲欧美| 久久久精品免费免费高清| 久久免费观看电影| 国产一区二区 视频在线| 欧美成人午夜精品| 久久久久精品人妻al黑| 国产成人一区二区三区免费视频网站| 国产亚洲精品一区二区www | 亚洲av欧美aⅴ国产| 老司机福利观看| 99国产精品一区二区三区| 国产精品久久久av美女十八| 亚洲av欧美aⅴ国产| 欧美国产精品一级二级三级| 精品国产超薄肉色丝袜足j| 人妻一区二区av| 国产精品偷伦视频观看了| 国产aⅴ精品一区二区三区波| 日韩中文字幕欧美一区二区| 91精品三级在线观看| 欧美成狂野欧美在线观看| 亚洲精品久久午夜乱码| 亚洲精品在线观看二区| 最近最新中文字幕大全免费视频| 高清在线国产一区| 美女视频免费永久观看网站| 一个人免费在线观看的高清视频| 亚洲成国产人片在线观看| 欧美乱码精品一区二区三区| 精品国产一区二区三区四区第35| 这个男人来自地球电影免费观看| 久久人人爽av亚洲精品天堂| 91字幕亚洲| 视频区欧美日本亚洲| 叶爱在线成人免费视频播放| 久久久久国产一级毛片高清牌| 精品一区二区三区av网在线观看 | 午夜福利在线免费观看网站| 中文字幕色久视频| 侵犯人妻中文字幕一二三四区| 亚洲第一青青草原| 午夜精品久久久久久毛片777| 久久青草综合色| 老熟妇仑乱视频hdxx| 亚洲精品国产精品久久久不卡| 免费在线观看黄色视频的| 色综合婷婷激情| 久久久国产欧美日韩av| 中亚洲国语对白在线视频| 一级毛片电影观看| 黄网站色视频无遮挡免费观看| 精品人妻熟女毛片av久久网站| 日韩有码中文字幕| 丰满少妇做爰视频| 国产一卡二卡三卡精品| 热99re8久久精品国产| 在线观看舔阴道视频| 一个人免费看片子| 欧美亚洲 丝袜 人妻 在线| 精品少妇久久久久久888优播| 三上悠亚av全集在线观看| 亚洲人成77777在线视频| 桃花免费在线播放| 少妇 在线观看| 宅男免费午夜| 日本vs欧美在线观看视频| 我要看黄色一级片免费的| 精品视频人人做人人爽| 99精品久久久久人妻精品| 大片免费播放器 马上看| 国产不卡一卡二| 免费人妻精品一区二区三区视频| 成人永久免费在线观看视频 | 精品第一国产精品| 人人妻,人人澡人人爽秒播| 少妇裸体淫交视频免费看高清 | 黄频高清免费视频| 久久热在线av| 亚洲美女黄片视频| 久久天躁狠狠躁夜夜2o2o| 淫妇啪啪啪对白视频| 黑人巨大精品欧美一区二区蜜桃| 国产免费现黄频在线看| 伦理电影免费视频| 老汉色∧v一级毛片| 精品人妻熟女毛片av久久网站| 亚洲伊人久久精品综合| 女同久久另类99精品国产91| 一二三四社区在线视频社区8| 久久精品国产亚洲av高清一级| 久久久久久久久久久久大奶| 久久精品熟女亚洲av麻豆精品| 欧美变态另类bdsm刘玥| 91成年电影在线观看| 欧美激情高清一区二区三区| 两性夫妻黄色片| 午夜日韩欧美国产| 人人妻,人人澡人人爽秒播| 人妻一区二区av| 十八禁网站免费在线| 国产亚洲精品久久久久5区| 一级片免费观看大全| 欧美亚洲 丝袜 人妻 在线| 中文字幕人妻熟女乱码| 777米奇影视久久| 一本久久精品| 中文亚洲av片在线观看爽 | 91成人精品电影| 成人亚洲精品一区在线观看| 视频区图区小说| 99九九在线精品视频| 日韩 欧美 亚洲 中文字幕| 高清黄色对白视频在线免费看| 国产男女内射视频| 午夜日韩欧美国产| 国产精品免费一区二区三区在线 | 性色av乱码一区二区三区2| 国产片内射在线| 少妇粗大呻吟视频| 久久人妻熟女aⅴ| 伦理电影免费视频| 午夜福利视频在线观看免费| 亚洲精品自拍成人| 一区二区三区国产精品乱码| 精品国产一区二区三区四区第35| 精品国产一区二区久久| 黑人巨大精品欧美一区二区mp4| 欧美日本中文国产一区发布| 欧美黄色淫秽网站| 国产一区二区三区视频了| 国产一区二区三区在线臀色熟女 | 国产成人av教育| 国产精品久久久久成人av| 他把我摸到了高潮在线观看 | 一区福利在线观看| 午夜老司机福利片| 一边摸一边做爽爽视频免费| 亚洲欧美一区二区三区黑人| 国产精品电影一区二区三区 | 男女之事视频高清在线观看| 看免费av毛片| 一区二区日韩欧美中文字幕| 十分钟在线观看高清视频www| 欧美成人午夜精品| 亚洲av日韩在线播放| 精品第一国产精品| avwww免费| 亚洲色图 男人天堂 中文字幕| 在线观看舔阴道视频| 国产欧美日韩一区二区三| 日本一区二区免费在线视频| 亚洲精品成人av观看孕妇| 性高湖久久久久久久久免费观看| 久久热在线av| 国产高清激情床上av| 亚洲人成电影免费在线| 又紧又爽又黄一区二区| 香蕉丝袜av| 成人18禁高潮啪啪吃奶动态图| 亚洲 国产 在线| 亚洲国产欧美一区二区综合| 在线观看舔阴道视频| 18禁国产床啪视频网站| av又黄又爽大尺度在线免费看| 天堂俺去俺来也www色官网| 亚洲性夜色夜夜综合| 精品第一国产精品| 日韩欧美一区视频在线观看| 99热网站在线观看| 亚洲精品中文字幕在线视频| 这个男人来自地球电影免费观看| 国产精品国产高清国产av | 日本一区二区免费在线视频| 亚洲精品在线美女| 国产黄色免费在线视频| 99国产精品一区二区三区| 在线天堂中文资源库| 亚洲欧洲日产国产| 国产xxxxx性猛交| 日本一区二区免费在线视频| 亚洲精品国产区一区二| a级毛片在线看网站| 99九九在线精品视频| 亚洲男人天堂网一区| 可以免费在线观看a视频的电影网站| 亚洲国产欧美日韩在线播放| 国产黄色免费在线视频| 欧美日韩av久久| 久久久国产精品麻豆| 丰满人妻熟妇乱又伦精品不卡| 日韩欧美国产一区二区入口| 亚洲一卡2卡3卡4卡5卡精品中文| 又紧又爽又黄一区二区| 国产精品电影一区二区三区 | 纵有疾风起免费观看全集完整版| 99国产精品99久久久久| 久久中文看片网| 亚洲国产欧美在线一区| 国产亚洲精品第一综合不卡| 99久久国产精品久久久| 黄网站色视频无遮挡免费观看| 日本a在线网址| 国产精品电影一区二区三区 | 日本黄色日本黄色录像| 国产成人啪精品午夜网站| 三上悠亚av全集在线观看| 亚洲一区二区三区欧美精品| 757午夜福利合集在线观看| 欧美日本中文国产一区发布| 制服人妻中文乱码| 久久精品aⅴ一区二区三区四区| 免费在线观看黄色视频的| 又紧又爽又黄一区二区| 午夜福利,免费看| 最新在线观看一区二区三区| 午夜精品国产一区二区电影| 日本av免费视频播放| 在线十欧美十亚洲十日本专区| 久久久国产精品麻豆| 精品一区二区三区av网在线观看 | 国产单亲对白刺激| 亚洲精品在线观看二区| 欧美人与性动交α欧美精品济南到| 老司机在亚洲福利影院| 国产福利在线免费观看视频| 咕卡用的链子| 91成年电影在线观看| 精品一区二区三区av网在线观看 | 亚洲第一欧美日韩一区二区三区 | 国产精品 国内视频| 欧美乱码精品一区二区三区| 黑丝袜美女国产一区| 正在播放国产对白刺激| 欧美变态另类bdsm刘玥| 色尼玛亚洲综合影院| 亚洲美女黄片视频| e午夜精品久久久久久久| 国产精品久久久久久精品古装| 满18在线观看网站| 老司机深夜福利视频在线观看| 亚洲中文av在线| 久久人人爽av亚洲精品天堂| 亚洲男人天堂网一区| 男女高潮啪啪啪动态图| 考比视频在线观看| www.精华液| 无遮挡黄片免费观看| 免费在线观看黄色视频的| 精品一区二区三区视频在线观看免费 | a级毛片黄视频| 免费观看人在逋| 精品国产国语对白av| 国产在线视频一区二区| 最近最新中文字幕大全电影3 | 日本黄色日本黄色录像| 免费在线观看日本一区| 91精品三级在线观看| 国产aⅴ精品一区二区三区波| 国产高清激情床上av| 欧美一级毛片孕妇| 欧美精品高潮呻吟av久久| 日本vs欧美在线观看视频| 午夜精品国产一区二区电影| bbb黄色大片| 另类精品久久| 成年人黄色毛片网站| 亚洲三区欧美一区| 日韩大码丰满熟妇| 一区二区三区激情视频| 18在线观看网站| 手机成人av网站| 中文字幕人妻丝袜制服| 一本综合久久免费| 99国产综合亚洲精品| netflix在线观看网站| 亚洲精品粉嫩美女一区| 日本av手机在线免费观看| 1024香蕉在线观看| 国产精品一区二区免费欧美| 免费在线观看完整版高清| 搡老乐熟女国产| 午夜福利乱码中文字幕| 亚洲欧洲日产国产| 少妇被粗大的猛进出69影院| 久久午夜综合久久蜜桃| 黄片小视频在线播放| 丰满人妻熟妇乱又伦精品不卡| 一区福利在线观看| 黄色丝袜av网址大全| 精品人妻在线不人妻| 亚洲综合色网址| 国产野战对白在线观看| 老司机午夜十八禁免费视频| 午夜免费成人在线视频| 国产精品亚洲av一区麻豆| 亚洲精品在线观看二区| 久久久久久亚洲精品国产蜜桃av| www.自偷自拍.com| 国产麻豆69| 黄色 视频免费看| 一区在线观看完整版| 国产精品久久久久久精品古装| 两性夫妻黄色片| 妹子高潮喷水视频| 久久九九热精品免费| 亚洲欧美激情在线| 国产黄频视频在线观看| 欧美 日韩 精品 国产| 一区在线观看完整版| 飞空精品影院首页| 午夜老司机福利片| 免费黄频网站在线观看国产| 中文字幕色久视频| h视频一区二区三区| 午夜成年电影在线免费观看| 欧美变态另类bdsm刘玥| 国产在线精品亚洲第一网站| www.999成人在线观看| 怎么达到女性高潮| 19禁男女啪啪无遮挡网站| 窝窝影院91人妻| av有码第一页| 精品国产国语对白av| 叶爱在线成人免费视频播放| 一本大道久久a久久精品| 国内毛片毛片毛片毛片毛片| 久久久久网色| 最新的欧美精品一区二区| 在线观看66精品国产| 十八禁人妻一区二区| 动漫黄色视频在线观看| 免费在线观看黄色视频的| 丝瓜视频免费看黄片| 黄色片一级片一级黄色片| 国产在线一区二区三区精| 精品久久久久久电影网| 搡老岳熟女国产| 国产男女内射视频| 亚洲成人国产一区在线观看| 亚洲 欧美一区二区三区| 黄色视频不卡| 亚洲av片天天在线观看| 精品免费久久久久久久清纯 | 两个人免费观看高清视频| √禁漫天堂资源中文www| 精品人妻在线不人妻| 成人影院久久| 99re6热这里在线精品视频| 高潮久久久久久久久久久不卡| 19禁男女啪啪无遮挡网站| 亚洲中文日韩欧美视频| 亚洲九九香蕉| 岛国在线观看网站| 高清毛片免费观看视频网站 | 免费一级毛片在线播放高清视频 | 侵犯人妻中文字幕一二三四区| 一区二区三区精品91| 999久久久精品免费观看国产| av福利片在线| 午夜福利视频在线观看免费| 午夜成年电影在线免费观看| 亚洲成av片中文字幕在线观看| 日日爽夜夜爽网站| 大陆偷拍与自拍| 免费日韩欧美在线观看| 啪啪无遮挡十八禁网站| 在线观看www视频免费| 亚洲人成电影免费在线| 无遮挡黄片免费观看| 美女扒开内裤让男人捅视频| 国产精品国产av在线观看| av天堂久久9| 老汉色av国产亚洲站长工具| 中文字幕高清在线视频| 肉色欧美久久久久久久蜜桃| 久久久久网色| 国产精品成人在线| 乱人伦中国视频| 日韩大片免费观看网站| 久久久精品免费免费高清| 久久久久久久久免费视频了| 露出奶头的视频| 一本久久精品| 麻豆av在线久日| 久久婷婷成人综合色麻豆| 亚洲国产av新网站| 天天影视国产精品| 伦理电影免费视频| 国产一区二区激情短视频| 午夜91福利影院| 伊人久久大香线蕉亚洲五| 热99re8久久精品国产| 精品亚洲成a人片在线观看| 美女扒开内裤让男人捅视频| 亚洲久久久国产精品| 亚洲国产欧美网| 黄频高清免费视频| 国产精品亚洲一级av第二区| 欧美乱妇无乱码| 亚洲欧美激情在线| 日韩大码丰满熟妇| 亚洲 国产 在线| 性少妇av在线| 极品少妇高潮喷水抽搐| 久久这里只有精品19| 搡老岳熟女国产| 91大片在线观看| 岛国毛片在线播放| 欧美人与性动交α欧美软件| 中国美女看黄片| 中文字幕制服av| 夜夜夜夜夜久久久久| 久久亚洲精品不卡| 日日夜夜操网爽| 免费女性裸体啪啪无遮挡网站| 久久国产精品影院| 人人妻人人澡人人看| 少妇精品久久久久久久| 在线播放国产精品三级| 在线天堂中文资源库| 中文字幕精品免费在线观看视频| 十八禁人妻一区二区| 成人国产一区最新在线观看| 国产精品免费大片| 国产欧美亚洲国产| 夜夜骑夜夜射夜夜干| 久久免费观看电影| 亚洲专区中文字幕在线| 欧美日韩亚洲国产一区二区在线观看 | 精品亚洲乱码少妇综合久久| 国产精品 欧美亚洲| 十八禁网站免费在线| 一本综合久久免费| 超色免费av| 亚洲专区中文字幕在线| 国产激情久久老熟女| 国产在线视频一区二区| 久久精品亚洲熟妇少妇任你| 成年人黄色毛片网站| 国产野战对白在线观看| 久久午夜综合久久蜜桃| 天天添夜夜摸| 精品少妇一区二区三区视频日本电影| 久久久精品国产亚洲av高清涩受| 久热这里只有精品99| 亚洲精品av麻豆狂野| 最近最新中文字幕大全免费视频| 美女视频免费永久观看网站| 婷婷丁香在线五月| 男女床上黄色一级片免费看| 亚洲自偷自拍图片 自拍| 这个男人来自地球电影免费观看| 大型黄色视频在线免费观看| 国产亚洲午夜精品一区二区久久| 女人久久www免费人成看片| 热99国产精品久久久久久7| 青草久久国产| videos熟女内射| 91大片在线观看| 大陆偷拍与自拍| 女同久久另类99精品国产91| av片东京热男人的天堂| 欧美激情 高清一区二区三区| 在线观看舔阴道视频| 又大又爽又粗| 国产亚洲精品第一综合不卡| 国产一区二区 视频在线| 久久中文字幕一级| 国产日韩欧美视频二区| 国产一区二区在线观看av| 欧美黑人精品巨大| 日日爽夜夜爽网站| 老鸭窝网址在线观看| 少妇粗大呻吟视频| 午夜久久久在线观看| 黄色视频不卡| 97人妻天天添夜夜摸| 免费在线观看黄色视频的| 成年动漫av网址| 国产精品电影一区二区三区 | 在线播放国产精品三级| 免费高清在线观看日韩| 亚洲熟女毛片儿| 日韩大片免费观看网站| 99国产精品免费福利视频| 老司机亚洲免费影院| 后天国语完整版免费观看| 亚洲国产毛片av蜜桃av| 久久久久久久久免费视频了| 亚洲精品中文字幕一二三四区 | 91麻豆av在线| 国产在线免费精品| 国产精品免费视频内射| 国产成人免费观看mmmm| 99精品欧美一区二区三区四区| 美国免费a级毛片| 一区二区三区精品91| 国产成人系列免费观看| 99精品在免费线老司机午夜| 电影成人av| 十八禁高潮呻吟视频| 久久热在线av| 满18在线观看网站| 亚洲专区中文字幕在线| 成人影院久久| 黄频高清免费视频| 午夜视频精品福利| 国产男女内射视频| 一个人免费在线观看的高清视频| 亚洲熟女毛片儿| av又黄又爽大尺度在线免费看| 19禁男女啪啪无遮挡网站| 国产精品98久久久久久宅男小说| 热99re8久久精品国产| 亚洲精品一二三| 国产在视频线精品| 久久午夜亚洲精品久久| 国产xxxxx性猛交| 美女福利国产在线| 亚洲精华国产精华精| videos熟女内射| 久久青草综合色| 69av精品久久久久久 | 国产成人啪精品午夜网站| 国产男女超爽视频在线观看| 妹子高潮喷水视频| 国产男女内射视频| 亚洲性夜色夜夜综合| 亚洲精品中文字幕在线视频| 欧美中文综合在线视频| 精品国产一区二区三区久久久樱花| 建设人人有责人人尽责人人享有的| 大码成人一级视频| 久久久国产精品麻豆| 国产亚洲一区二区精品| 亚洲avbb在线观看| 人人妻人人澡人人爽人人夜夜| 日韩一区二区三区影片| 久久久久视频综合| 成人免费观看视频高清| 精品国产亚洲在线| 91麻豆av在线| 露出奶头的视频| 999久久久精品免费观看国产| 国产激情久久老熟女| 叶爱在线成人免费视频播放| 少妇粗大呻吟视频| 亚洲欧洲精品一区二区精品久久久| 久久人妻熟女aⅴ|