馮艷霞，張志紅，張少強(qiáng)

(天津師范大學(xué)計(jì)算機(jī)與信息工程學(xué)院，天津300387)

(*通信作者電子郵箱zhangshaoqiang@tjnu.edu.cn)

0 引言

基于新一代基因測(cè)序的高通量技術(shù)，特別是染色質(zhì)免疫共沉淀的高通量測(cè)序(ChromatinImmunoprecipitation Sequencing，ChIP-Seq)技術(shù)［1］，在全基因組范圍內(nèi)的使用，很大程度上改變了生物學(xué)家定位大規(guī)模真核基因組順式調(diào)控元件和模體的方式［2］。順式調(diào)控模體是同一轉(zhuǎn)錄因子結(jié)合位點(diǎn)的序列集合。模體發(fā)現(xiàn)問題一直是生物信息學(xué)比較重要的研究課題之一。模體所蘊(yùn)含的遺傳信號(hào)，對(duì)于在全基因組范圍內(nèi)研究基因的轉(zhuǎn)錄與調(diào)控機(jī)制具有重大意義。在ChIP-Seq技術(shù)之前的模體發(fā)現(xiàn)，主要是通過比較同源基因或共調(diào)控基因非編碼區(qū)域的保守性來預(yù)測(cè)調(diào)控元件模體，其主要方法包括:通過多序列比對(duì)尋找保守序列片段;通過位置賦權(quán)矩陣更新來遍歷調(diào)控元件;運(yùn)用統(tǒng)計(jì)取樣的估計(jì)方法等。這些方法均取得了較好的預(yù)測(cè)結(jié)果，并應(yīng)運(yùn)而生了大量模體發(fā)現(xiàn)算法和工具，例如多重期望最大模體提取(Multiple Expectation maximization for Motif Elicitation，MEME)算法［3］、MotifSampler算法［4］、Weeder算法［5］等，為生物學(xué)家對(duì)真核生物全基因組［6］的解讀和研究提供了有利保障，推動(dòng)著生物信息學(xué)的蓬勃發(fā)展。模體發(fā)現(xiàn)算法的工作流程如圖1所示。

由于ChIP-Seq單次實(shí)驗(yàn)產(chǎn)生了海量的堿基讀序(reads)，通過序列比對(duì)拼接及結(jié)合峰預(yù)測(cè)［7］后仍有成千上萬條序列，這比傳統(tǒng)的同源基因或共調(diào)控基因序列規(guī)模大很多。為了處理如此大規(guī)模的數(shù)據(jù)，近來的算法主要采用“有區(qū)別的”(discriminative)模型［8］策略，主要包括:一是選擇數(shù)據(jù)集的子集作為實(shí)驗(yàn)的輸入;二是在消耗較大數(shù)據(jù)存儲(chǔ)空間的情況下，運(yùn)用新型數(shù)據(jù)結(jié)構(gòu)來節(jié)省運(yùn)行時(shí)間;三是將實(shí)驗(yàn)數(shù)據(jù)集分割成兩個(gè)不等的部分，即實(shí)驗(yàn)組和對(duì)照組，僅將實(shí)驗(yàn)組作為模體發(fā)現(xiàn)算法的核心運(yùn)行數(shù)據(jù)，然后根據(jù)實(shí)驗(yàn)組發(fā)現(xiàn)的模體，與對(duì)照組進(jìn)行比對(duì)。上述第一種方法數(shù)據(jù)損失較多，要預(yù)測(cè)的結(jié)合位點(diǎn)模體信息在子集中可能不保守，從而大大影響算法準(zhǔn)確性，例如獨(dú)立成分分析(Independent Component Analysis，ICA)和套層抽樣相結(jié)合的的模體發(fā)現(xiàn)算法——NestedMICA(Nested saMpling and ICA)［9］;第二種方法的主要問題是隨著數(shù)據(jù)規(guī)模增加，會(huì)成倍增加內(nèi)存需求和時(shí)間復(fù)雜性，從而降低了算法運(yùn)行速度，例如有區(qū)別的正則表達(dá)式模體提取(Discriminative Regular Expression Motif Elicitation，DREME)算法［10］;相比較而言，第三種方法更注重序列數(shù)據(jù)集的分隔方法，但不恰當(dāng)?shù)臄?shù)據(jù)分隔可能導(dǎo)致錯(cuò)誤的結(jié)果，例如Amadeus［11］、Trawler［12］。特別是，隨著 ChIP-Seq 數(shù)據(jù)的海量增長，其噪聲也越來越多，從而導(dǎo)致預(yù)測(cè)結(jié)果具有較高的錯(cuò)誤發(fā)生率(False Discovery Rate，F(xiàn)DR)。此外，模體發(fā)現(xiàn)算法在追求時(shí)間與空間的高效外，還應(yīng)準(zhǔn)確識(shí)別出模體的真實(shí)長度。

為了解決上述問題，本文設(shè)計(jì)了一個(gè)基于ChIP-Seq數(shù)據(jù)集的順式調(diào)控模體發(fā)現(xiàn)算法——FisherNet，該算法既能提高模體預(yù)測(cè)的準(zhǔn)確率，又能提高算法運(yùn)行速度。

圖1 模體發(fā)現(xiàn)流程Fig.1 Flow chart of motif finding

1 數(shù)據(jù)來源與預(yù)處理

為了檢驗(yàn)算法的優(yōu)越性，F(xiàn)isherNet算法使用與DREME［10］算法相同的數(shù)據(jù)集，包括13個(gè)小鼠胚胎干細(xì)胞(Mouse Embryonic Stem Cell，mESC)ChIP-Seq 數(shù)據(jù)集［13］、3 個(gè)小鼠紅細(xì)胞數(shù)據(jù)集和1個(gè)人類淋巴母細(xì)胞系數(shù)據(jù)集。mESC數(shù)據(jù)集為12個(gè)轉(zhuǎn)錄因子，它是維持多能性和CTCF的關(guān)鍵。小鼠紅細(xì)胞的數(shù)據(jù)集為 Gata1［14］、Klf1［15］和 Tal1，是紅細(xì)胞生成的關(guān)鍵調(diào)節(jié)者。人體淋巴細(xì)胞數(shù)據(jù)集是維生素D受體包括在使用骨化三醇的刺激前和之后的ChIP-Seq數(shù)據(jù)。另外，基于人的 ENCODE(Encyclopedia of DNA Elements)［16］數(shù)據(jù)庫的ChIP-Seq數(shù)據(jù)集規(guī)模超大，本文選取其中的10個(gè)有已知模體的轉(zhuǎn)錄因子(ATF1、ATF3、CTCF、E2F4、EBF1、ECF1、POL2、ZNF274、USF2、YY1)數(shù)據(jù)集進(jìn)行比較。這些數(shù)據(jù)集源文件格式為FASTQ格式。

2 算法的主要步驟

本文算法的流程如圖2所示。

2.1 生成背景數(shù)據(jù)集

FisherNet算法采用了三階馬爾可夫鏈［17］生成背景數(shù)據(jù)，即遍歷實(shí)驗(yàn)組ChIP-Seq數(shù)據(jù)集統(tǒng)計(jì)每三個(gè)堿基短序列(共43=64種)后出現(xiàn)各個(gè)堿基的頻率得到概率轉(zhuǎn)移矩陣。然后用該轉(zhuǎn)移矩陣生成與實(shí)驗(yàn)組數(shù)據(jù)中序列長度和數(shù)目相同的背景序列數(shù)據(jù)集。

2.2 遍歷短序生成哈希表

為了統(tǒng)計(jì)每個(gè)長度為k的短序(k-mer)［18］在各個(gè)序列是否存在，F(xiàn)isherNet分別遍歷實(shí)驗(yàn)數(shù)據(jù)集和背景數(shù)據(jù)集中的每條序列并用兩個(gè)哈希表分別存儲(chǔ)所有的k-mer信息。哈希表的鍵為k-mer，鍵值為該k-mer所在序列編號(hào)集合。因?yàn)檎婧松锏慕Y(jié)合位點(diǎn)長度主要介于4～8個(gè)堿基長。FisherNet算法默認(rèn)k從4到8個(gè)堿基為窗口長度分別進(jìn)行遍歷。在遍歷過程中，同時(shí)考慮k-mer在序列的互補(bǔ)鏈?zhǔn)欠翊嬖?，若存在，則將該序列的編號(hào)也放到該k-mer鍵值對(duì)應(yīng)的集合中。遍歷過程如圖3所示。

圖2 本文算法流程Fig.2 Flow chart of proposed algorithm

2.3 根據(jù)費(fèi)舍爾精確檢驗(yàn)計(jì)算P值

FisherNet算法采用費(fèi)舍爾精確檢驗(yàn)分析每個(gè)k-mer的富集情況［19］。對(duì)于每個(gè)k-mer，通過哈希表查詢記錄下它在實(shí)驗(yàn)數(shù)據(jù)集和背景數(shù)據(jù)集中出現(xiàn)的序列條數(shù)。如表1所示，a代表某k-mer在實(shí)驗(yàn)組數(shù)據(jù)集中出現(xiàn)過的序列數(shù)，b代表該短序在背景數(shù)據(jù)集中出現(xiàn)過的序列數(shù)，相對(duì)地，c和d分別代表該短序在兩組實(shí)驗(yàn)數(shù)據(jù)集中未出現(xiàn)的序列數(shù)。

根據(jù)表1通過超幾何分布公式計(jì)算P值［20］為:

若P值大于意義臨界值0.05(假設(shè)檢驗(yàn)常規(guī)設(shè)P值的閾值為0.05，在算法中該閾值作為輸入?yún)?shù)可由使用者自主調(diào)整)，說明該條DNA短序列存在基因信息與背景差異較小，則對(duì)其不再考慮。上述計(jì)算P值的式(1)，對(duì)于小規(guī)模數(shù)據(jù)集效果良好，對(duì)于大型ChIP-Seq數(shù)據(jù)集，則階乘n!會(huì)是超大的數(shù)，甚至超出計(jì)算機(jī)整數(shù)計(jì)數(shù)范圍。對(duì)此，F(xiàn)isherNet引入了斯特林(Stirling)公式，對(duì)組合數(shù)的計(jì)算先取其對(duì)數(shù)ln，再取指數(shù)。斯特林如式(2)所示，與ln C推導(dǎo)如式(3)、(4)所示:

則P值公式就等價(jià)于:

對(duì)式(1)兩邊取對(duì)數(shù)得:

對(duì)于每個(gè)k-mer，用式(6)計(jì)算其P值，若P值小于0.05，以此k-mer為鍵，其P值為鍵值保存到一個(gè)新的哈希表。通過比較各k-mer的鍵值，對(duì)k-mer進(jìn)行排序得到候選k-mer集合。

表1 k-mer(模體)序列在實(shí)驗(yàn)組和背景組的數(shù)據(jù)統(tǒng)計(jì)Tab.1 Statistics of k-mer(motif)sequence in experimental group and background group

2.4 構(gòu)建初始模體

在候選k-mer集合中選取P值最小的k-mer作為種子，依次向下遍歷與種子長度相同且只有一個(gè)位置不同的k-mer，將其合并形成初始模體。每合并一個(gè)k-mer就要通過計(jì)算出現(xiàn)模體和未出現(xiàn)模體的序列數(shù)來更新費(fèi)舍爾精確檢驗(yàn)表并計(jì)算合并后模體的P值，若P值 ＞0.05，則停止k-mer合并。統(tǒng)計(jì)初始序列上各堿基的出現(xiàn)頻率得到位置頻率矩陣P=(P(b，i))4×k［21］，其中 P(b，i)表示堿基 b 在位置 i出現(xiàn)的頻率;然后再根據(jù)式(7)計(jì)算出對(duì)應(yīng)的位置賦權(quán)矩陣W =(W(b，i))4×k［22］，其中:P(b)是堿基 b 在實(shí)驗(yàn)數(shù)據(jù)集上出現(xiàn)的頻率(為了能夠計(jì)算對(duì)數(shù)值，當(dāng)P(b，i)=0，令P(b，i)=10－5)。構(gòu)建初始模體的流程如圖4中第① ～③步所示。

2.5 生成最終模體

用上述生成的位置賦權(quán)矩陣掃描候選k-mer集合中的每一個(gè)相同長度的k-mer={b1，b2，…，bk}，為了衡量一個(gè)k-mer是否屬于模體，本文用式(8)(公式中的Smax、Smin分別為位置賦權(quán)矩陣中每列元素的最大值的和與最小值的和)計(jì)算k-mer的相對(duì)分?jǐn)?shù)值R并只保留R≥0.80(根據(jù)已知模體統(tǒng)計(jì)有80% 的位置具有保守性)的k-mer，以此生成最終模體。如圖4中第④到⑤步所示。

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

FisherNet算法運(yùn)用mESC的ChIP-Seq數(shù)據(jù)集和一個(gè)意義臨界值對(duì)結(jié)果進(jìn)行具體分析。如表2為FisherNet算法預(yù)測(cè)出的模體發(fā)現(xiàn)情況。將算法發(fā)現(xiàn)的模體序列與Tomtom數(shù)據(jù)庫［23］中的真實(shí)模體進(jìn)行比較，針對(duì)Tomtom標(biāo)準(zhǔn)數(shù)據(jù)庫中已有的13個(gè)重要模體，F(xiàn)isherNet算法發(fā)現(xiàn)了其中的10個(gè)。

此外，針對(duì)mESC ChIP-Seq數(shù)據(jù)集，本文分別運(yùn)行了其他六種被廣泛應(yīng)用且能發(fā)現(xiàn)輔調(diào)控因子模體的模體發(fā)現(xiàn)算法，即 MEME［3］、Weeder［5］、NestedMICA［9］、DREME［10］、Amadeus［11］、Trawler［12］，并將這些算法的預(yù)測(cè)結(jié)果與 FisherNet順式調(diào)控模體算法相比較，統(tǒng)計(jì)結(jié)果如表3所示。其中:N為各算法發(fā)現(xiàn)的模體總數(shù);S為與Tomtom標(biāo)準(zhǔn)數(shù)據(jù)庫中相匹配的模體數(shù);C為算法發(fā)現(xiàn)的輔調(diào)控因子模體數(shù)。在發(fā)現(xiàn)轉(zhuǎn)錄因子模體上FisherNet與 DREME、Weeder、NestedMICA 和 MEME 均達(dá)到最高，但在發(fā)現(xiàn)輔調(diào)控因子模體上FisherNet算法要優(yōu)于其他算法。這說明該算法在保證發(fā)現(xiàn)模體的精確度的同時(shí)，能發(fā)現(xiàn)更多的輔調(diào)控因子模體，F(xiàn)isherNet總體能夠達(dá)到80%(即(6.1+10)/20)的精確度。表3最后一列為各算法對(duì)13個(gè)mESC數(shù)據(jù)集的平均運(yùn)行時(shí)間，可以看出，F(xiàn)isherNet是所有比對(duì)算法中平均運(yùn)行時(shí)間最短的，而NestedMICA和Weeder是運(yùn)行時(shí)間最長的兩個(gè)算法。

表2 小鼠的ES細(xì)胞ChIP-Seq數(shù)據(jù)集發(fā)現(xiàn)的模體Tab.2 Found motif in mouse ES cells ChIP-Seq dataset

表3 不同模體發(fā)現(xiàn)算法結(jié)果統(tǒng)計(jì)Tab.3 Results statistics of different motif finding algorithms

由于DREME是MEME開發(fā)者設(shè)計(jì)的改進(jìn)算法，為了進(jìn)一步檢驗(yàn)除了MEME、Weeder和NestedMICA之外的四種算法的效率，分別取序列規(guī)模為1兆、2兆、3兆和4兆堿基對(duì)的數(shù)據(jù)集進(jìn)行運(yùn)算，圖5為四種算法的運(yùn)行時(shí)間評(píng)估折線圖。由圖5可以看出，與其他算法相比，隨著數(shù)據(jù)集序列堿基規(guī)模的增加，F(xiàn)isherNet算法運(yùn)行時(shí)間最短且增加幅度不大，充分說明了該算法在運(yùn)行效率上的優(yōu)越性。將上述算法應(yīng)用到選取的ENCODE［16］數(shù)據(jù)集中，發(fā)現(xiàn)每一個(gè)真實(shí)模體均存在于FisherNet算法輸出的前10個(gè)模體中，而其余算法的前10個(gè)模體中含有真實(shí)模體的準(zhǔn)確度均小于80%。

圖5 不同算法運(yùn)行時(shí)間比較Fig.5 Running time comparison of different motif finding algorithms

4 結(jié)語

本文基于ChIP-Seq數(shù)據(jù)的模體發(fā)現(xiàn)，以ChIP-Seq結(jié)合峰數(shù)據(jù)集為算法輸入數(shù)據(jù)，提出一個(gè)全新的順式調(diào)控模體發(fā)現(xiàn)算法(本文的FisherNet由 Perl語言編寫，代碼下載地址為http://www.escience.cn/people/sqzhang/program.html.)。該算法運(yùn)用費(fèi)舍爾假設(shè)檢驗(yàn)來尋找區(qū)分度高的k-mer并通過位置賦權(quán)矩陣進(jìn)行優(yōu)化篩選。與其他多種模體發(fā)現(xiàn)算法進(jìn)行比較，驗(yàn)證了該算法大大提高了模體發(fā)現(xiàn)的數(shù)目和準(zhǔn)確度，成功預(yù)測(cè)到許多核心模體及輔調(diào)控因子模體。該算法目前只對(duì)幾種代表性的ChIP-Seq數(shù)據(jù)進(jìn)行了驗(yàn)證，對(duì)其他物種的數(shù)據(jù)是否具有魯棒性有待進(jìn)一步驗(yàn)證。