基于顯露子串挖掘的基因序列模體識別算法

張懿璞+閆茂德++侯俊++闞丹會

摘 要： 染色質(zhì)免疫共沉淀技術(shù)將模體識別問題拓展到了全基因組范圍，但因數(shù)據(jù)量過大，傳統(tǒng)的模體識別算法往往運算過慢從而無法很好地解決此問題。為了解決傳統(tǒng)算法的缺點，提出一種用于ChIP?seq數(shù)據(jù)的替換顯露子串尋找問題的算法FastESE，通過測試集和控制集的比對找出顯露子串并搜索其（l，d）替換實例組成相應(yīng)的位置概率矩陣，再使用權(quán)重信息量對這些子串進行聚類，最終找出集合中的替換顯露子串。使用真實的ChIP?seq數(shù)據(jù)對該研究算法進行有效性驗證，實驗結(jié)果表明，F(xiàn)astESE可以在合理時間內(nèi)有效解決ChIP?seq中的模體識別問題。

關(guān)鍵詞： 染色質(zhì)免疫共沉淀； 顯露子串； 模體識別； FastESE

中圖分類號： TN911?34； TP301.6 文獻標(biāo)識碼： A 文章編號： 1004?373X（2017）12?0006?05

Abstract： Recently， the development of chromatin immunoprecipitation technique has extended the motif identification problem to the genome?wide range， but the traditional motif identification algorithms runs too slowly and hard to solve this large?scale data problem. In order to solve the shortcomings of the traditional algorithms， a substituted emerging substring search algorithm named FastESE applied to ChIP?seq data is proposed in this research. The emerging substrings are found out by comparing the test dataset and the control dataset， and then its substituted instances are searched to constitute the corresponding position probabilistic matrix. The weighted information content is adopted to cluster these substrings， and Finally， discover the substituted emerging substrings. The effectiveness of proposed algorithm was verified with the real ChIP?seq data. The experimental results show that the FastESE can deal with the motif identification problem in the ChIP?seq data in a proper time.

Keywords： chromatin immunoprecipitation； emerging substring； motif identification； FastESE

0 引 言

模體識別問題是研究基因序列的一個重要并且有挑戰(zhàn)性的問題，一直以來受到計算機學(xué)家和生物學(xué)家的重點研究。隨著染色質(zhì)免疫共沉淀技術(shù)（Chromatin Immunoprecipitation，ChIP）的發(fā)展[1?2]，染色質(zhì)免疫共沉淀技術(shù)與高通量測序方法結(jié)合的ChIP?seq技術(shù)將模體識別問題拓展到了全基因組范圍。而這類ChIP?seq數(shù)據(jù)量比起傳統(tǒng)的模體識別數(shù)據(jù)，數(shù)據(jù)規(guī)模大大增加[3]，這也就使得傳統(tǒng)的模體識別算法難以很好地解決ChIP?seq數(shù)據(jù)的模體識別問題[4]。

近些年，一些傳統(tǒng)算法也提出了用于ChIP數(shù)據(jù)的版本，如HMS，STEME，Weeder[5?7]等。這些算法可以解決部分ChIP?seq模體識別問題，但仍存在不足，例如HMS在訓(xùn)練時只能使用部分得分較高的序列；而STEME使用后綴樹加速了期望最大化求精，但只能用于尋找較短的模體；Weeder只能將單個模體作為每次運行的輸出結(jié)果，不符合實際需求。

針對ChIP?seq數(shù)據(jù)含有序列條數(shù)較多（上千條），但長度相對較短的特征（通常為幾百個堿基）。一些研究者將模體識別問題轉(zhuǎn)化為替換顯露子串挖掘問題[8]（Substituted emerging substrings mining problem），即從背景成分較少，相對較為干凈的大量序列中，區(qū)分出一些具有顯著性特征的子串。

本文中，針對上述ChIP數(shù)據(jù)的特點和傳統(tǒng)模體識別算法存在的缺陷，提出一種新的算法FastSES（Fast Substituted Emerging Substrings finding algorithm），用于ChIP?seq數(shù)據(jù)中替換顯露子串的搜索，F(xiàn)astSES通過含有顯露子串的測試集和不含顯露子串的控制集進行比對，篩選出頻率較高的顯露子串并生成對應(yīng)的概率分布矩陣，再對其進行聚類尋找最終的替換顯露子串。使用老鼠胚胎干細(xì)胞上的多組ChIP?seq數(shù)據(jù)對本文算法進行測試，實驗結(jié)果表明本文算法可以在合理時間內(nèi)有效找出ChIP?seq序列中的模體。

1 相關(guān)工作

1.1 問題定義

令Σ={k1，k2，…，km}為序列生成所使用的字母表，由Σ生成的序列集合S={S1，S2，…，St}，集合S中任一序列Si=，其長度均為n（s∈Σ）。已知長度為l的子串A=與長度為n的序列S=，如果對于，都有=，那么A就為S的一個子串，記為。如果A≠S，記為。

（1） 顯露子串挖掘問題（Emerging substrings mining problem）。已知由字符集Σ生成的測試序列集合St和控制序列集合Sc，令λf為頻率閾值，λg （λg>1）為增長率閾值。則對于一個輸入集合中的子串u，當(dāng)f（u，St）≥λf并且g（u，St，Sc）≥λg時，稱該字符串為一個顯露子串。這里f（u，S）表示子串u在集合S中發(fā)生的頻率，表示子串u在測試集合St中相對于控制集合Sc的增長率。g（u，St，Sc）的值越大，子串u在兩個集合中的相對分辨度就越高。顯露子串挖掘問題就是通過搜索序列集合中的每一個子串，尋找那些同時滿足f（u，St）≥λf和g（u，St，Sc）≥λg的子串。

（2） 替換顯露子串挖掘問題（Substituted emerging substrings mining problem）。對一個由字符集Σ生成的長度為l的字符串u，如果存在與其長度相同的子串u′，u′與u在至多d（d

1.2 數(shù)學(xué)模型

顯露子串通?？梢杂胒wk表示字符k出現(xiàn)在子串第w個位置的頻率，f0k為背景成分中出現(xiàn)字符k的概率。對于序列中任一長度為l的子串A=，第w個位置字符aw的對數(shù)似然值即為：

為了更好地表示顯露子串，在傳統(tǒng)的模體成分和背景成分之外還考慮到各位置間的內(nèi)在依賴關(guān)系，也就是某些連續(xù)字符組合的頻率脫離了獨立模體分布的期望概率[9]。如在連續(xù)位置出現(xiàn)的字符“AC”的概率與第一個位置出現(xiàn)字符“A”，第二個位置出現(xiàn)字符“C”的聯(lián)合概率有明顯的差異，即認(rèn)為這兩個位置存在內(nèi)在依賴關(guān)系。令Φw，w+1表示子串內(nèi)第w和第w+1個位置出現(xiàn)字符aw，aw+1的內(nèi)在依賴概率，其對數(shù)似然值為：

式中，Φ0表示背景概率。由式（2）可計算子串x的對數(shù)似然比（Log?Likelihood Ratio，LLR）：

式（4）為x各位置字符獨立時的概率，而式（5）為x各位置字符考慮內(nèi)在依賴關(guān)系時的概率。

2 算法描述

對于ChIP?seq數(shù)據(jù)中的替換顯露子串挖掘問題，采用差異模體發(fā)現(xiàn)方法，也就是輸入含有替換顯露子串的測試集和不含替換顯露子串的控制集兩個集合，首先找到顯露子串問題，再對替換子串進行拓展聚類，尋找相應(yīng)的替換顯露子串。通常，在ChIP?seq數(shù)據(jù)中，測試集是ChIP?seq的峰值區(qū)域，控制集可以是非峰值區(qū)域或者不同ChIP?seq數(shù)據(jù)中的類似序列?；谏鲜龇椒?，F(xiàn)astESE算法由詞頻分析、位置頻率矩陣構(gòu)建、聚類三個階段構(gòu)成。

2.1 詞頻統(tǒng)計

對于長度為l的子串（l?mer），其在數(shù)據(jù)集中可能出現(xiàn)的種類有4l個，那么首先在4l個子串中搜索顯露子串，即在測試集中頻繁出現(xiàn)而在控制集中很少出現(xiàn)的那些子串。

令S1表示輸入的測試集，通常為含有模體的ChIP數(shù)據(jù)；S2為控制集，通常使用不含模體的背景序列作為參照。由于顯露子串長度l未知，一般的，選取長度范圍為6～12的子串，統(tǒng)計其相應(yīng)數(shù)量并將每一個固定長度的子串存儲在O（4l）的空間中。為了方便理解，在表1中舉例表示對每一個l?mer計數(shù)并搜索顯露子串。對于同時滿足f（u，St）>λf和g（u，St，Sc）>λg的子串，如“TCTGAG”，即為顯露子串。

如果輸入的測試集和控制集均含有長度為n的t條序列，那么進行詞頻統(tǒng)計尋找顯露子串的計算復(fù)雜度僅為O（ntl），這僅與傳統(tǒng)訓(xùn)練方法中一次迭代的計算量相同。通過本步驟，可以快速地移除背景成分的干擾，大大降低運算量。

表1 顯露子串搜索樣例

2.2 構(gòu)建位置頻率矩陣

由于顯露子串只是對子串進行的基于字符的一致序列描述，并無法反映其統(tǒng)計顯著性。需要生成其相應(yīng)的位置頻率矩陣，從而進一步地衡量其統(tǒng)計顯著性。在本步驟中，通過搜索每一個顯露子串的替換（l，d）實例并計算每個替換實例的z值進行篩選，從而生成對應(yīng)的位置頻率矩陣。這里使用的z值是基于超幾何概率分布（Hyper?geometric Probability Distribution）的估計得到[10]：

式中：A1和A2分別表示某一個子串在集合S1和S2中出現(xiàn)的數(shù)量；C1和C2分別表示集合S1和S2中子串的總數(shù)量；。同樣通過表2來舉例說明位置頻率矩陣生成的過程。

在表2中，顯露子串“ACCACGTG”在S1和S2中分別出現(xiàn)119次和9次，已知l=8，選取d=1，即（l，d）=（8，1），在24個替換子串中篩選顯露子串“ACCACGTG”的顯露替換子串。接下來，分別計算24個替換子串的z值，并選擇z>1.643的子串為合格的顯露替換子串（如替換子串的z值都小于1.643，取z值為正即可）。將合格的顯露替換子串與顯露子串在測試集和控制集中每個位置上出現(xiàn)字符的數(shù)量相加，即可得到顯露子串在測試集和控制集中對應(yīng)的位置計數(shù)矩陣N1，N2。注意到，由于測試集和控制集的大小不同，那么控制集中子串的數(shù)量需要使用進行調(diào)整；并且，比起初始的搜索空間4l = 65 536，使用（l，d）進行替換子串的搜索可以使得計算量大大降低。

考慮到測試集中包含有所要尋找的顯露子串和背景序列，而控制集是完全由背景序列構(gòu)成，那么為了最大限度地消除測試集中顯露子串受到背景序列的干擾，使用控制集對測試集的預(yù)測結(jié)果進行糾正。也就是，將N2看作是N1中背景序列的期望計數(shù)矩陣，則所要尋找的顯露子串的位置計數(shù)矩陣N=max（N1-N2，0）。在N中的每個位置，還加入1%的偽計數(shù)防止零概率的出現(xiàn)，通過N即可得到顯露子串的概率分布與字符的內(nèi)在依賴關(guān)系分布。

2.3 聚 類

由各顯露子串的概率分布，就可以通過計算對數(shù)似然比（LLR）來評價其顯著性，見式（3）。而針對那些較為相近的顯露子串，使用基于LLC的WIC值（Weighted Information Content）來度量子串間的相似度[11]。WIC由子串α和子串β定義為：

（7）

式中：δ表示模體α的第δ個位置；ψ為模體β的第ψ個位置；ξ一般取值[11]為2.5；LLC（αδ）表示模體α第δ個位置的對數(shù)似然值；DIC（Differential Information Content）表示信息量差異，定義如下：

（8）

當(dāng)一個新的顯露子串被找到時，首先根據(jù)式（7）計算WIC值來檢查是否有已預(yù)測出的相似顯露子串，如果有，比較其LLR，保留LLR高的；如果沒有，用新顯露子串代替原顯露子串中LLR最小的，如此進行直至沒有新顯露子串被發(fā)現(xiàn)。FastESE算法流程如下：

輸入：測試集S1控制集S2

輸出：模體集C

For l ← lmin to lmax do

For each l?mer of 4l l?mers： u do

if f（u，St）≥λf && g（u，St，Sc）≥λg then

add u to Xs

For each l?mer of Xs：x do

For each d ← 1 to dmax do

calculate z?score of each neighborhood instance x′

if z（x′）>1.643 then

Add x′ to Xq

use x and Xq to construct Θ and Φ

add Θ to set A and add Φ to set B

For each Θ of A do

use Θ and corresponding Φ to compute LLC and WIC.

add xmotif formed by Θ of top LLR score to C

return C

其中，第1～4行為詞頻統(tǒng)計，尋找顯露子串的過程；第5～11行為通過顯露子串，生成其相應(yīng)位置頻率矩陣的過程；第12～15行為聚類并輸出最終結(jié)果。

3 實驗分析

為了測試本文算法性能，在此選擇了12組具有物種多樣性特征的老鼠胚胎干細(xì)胞（mES cells）ChIP?seq數(shù)據(jù)[12]，具體信息如表3所示。

其中Nanog，Oct4，Sox2，Esrrb和Zfx為自我復(fù)制的調(diào)控區(qū)域提取出的數(shù)據(jù)；Tcfcp2l1是胚胎干細(xì)胞中的優(yōu)先調(diào)控區(qū)域提出的數(shù)據(jù)；Klf4，cMyc，nMyc是重組因子數(shù)據(jù)；Smad1和STAT3對信號通路有著重要意義；而CTCF是轉(zhuǎn)錄隔離的重要成分。通過這些老鼠胚胎干細(xì)胞的ChIP?seq數(shù)據(jù)，選取峰值為中心，左右各100bp的序列片段為測試集，而在測試集序列兩端向外再延伸400bp后，選取長度500bp的序列片段作為控制集合。

FastESE算法運行參數(shù)設(shè)置如下：子串搜索長度范圍為6～12，當(dāng)每條序列中含有一個替換顯露子串時λf=0.8Pocc；當(dāng)每條序列中含有一個或不含替換顯露子串時λf=0.6Pocc；當(dāng)每條序列中含有零個或多個替換顯露子串時λf=1.2Pocc。Pocc=0.2，0.5，0.8分別表示高保守性、中保守性和低保守性，λg=2。搜索替換顯露子串的（l，d）包括（6，1），（7，1），（8，1），（9，2），（10，2），（11，2）和（12，3）。z=1.643為默認(rèn)閾值。使用Matlab語言在Windows系統(tǒng)環(huán)境下實現(xiàn)FastESE算法，測試運行環(huán)境為2.67 Hz CPU和4 GB內(nèi)存。

為了更好地顯示算法性能，使用Chen等人使用的WEEDER算法所公布的模體與本文方法進行比較。為了掌握預(yù)測模體和所公布數(shù)據(jù)間的相似性[13]，采用將統(tǒng)計顯著性表現(xiàn)為LOGO圖的方式，這也是現(xiàn)在被廣為使用的一種對比方式。圖1為算法發(fā)現(xiàn)的主要模體和Chen公布模體的對比，可以看到本文的算法可以有效地尋找到這些真實數(shù)據(jù)中的主要模體，LOGO圖顯示了較好的準(zhǔn)確率。

FastESE預(yù)測的模體LOGO圖

除準(zhǔn)確性外，還與現(xiàn)在流行的模體識別算法在時間效率上進行了比較，見圖2。

從圖2的結(jié)果可以看到，本文算法比起Weeder，HMS和ChIPMunk等流行算法，具有更好的效率，在處理更大規(guī)模的基因序列時所用時間最短。與FastESE相比，Weeder算法在設(shè)計時只是針對傳統(tǒng)模體識別問題，所以在處理上千條序列的數(shù)據(jù)時效率很低；而HMS在處理ChIP數(shù)據(jù)時只是選取部分序列進行模體識別，不能完全處理全部的數(shù)據(jù)，這樣就導(dǎo)致處理大規(guī)模序列時訓(xùn)練時間較長、迭代緩慢，從而影響效率；ChIPMunk同樣也存在收斂緩慢、訓(xùn)練時間較長的缺點。而本文算法在處理上百兆數(shù)據(jù)時，僅需要幾分鐘，如在處理Sox2數(shù)據(jù)時，如果?。╨，d）=（8，1），本文算法需要536 s，而（l，d）=（10，1）時，需運行620 s。

此外，將Nanog與Sox2或Oct4數(shù)據(jù)集對比，還使用FastESE算法找出了部分顯著性較弱的復(fù)合模體。將Nanog的數(shù)據(jù)作為測試集，Sox2或Oct4數(shù)據(jù)作為控制集作為輸入，Nanog VS Oct4的結(jié)果為“GCAATCAA”，Nanog VS Sox2的結(jié)果則為“CCATCAA”。也就是說，在這個相似區(qū)域中，Nanog模體與Sox2和Oct4轉(zhuǎn)錄因子中的一個或兩個間接綁定。這些顯著性較弱的復(fù)合模體與Chen等人公布的Sox2和Oct4數(shù)據(jù)中有相似的結(jié)合區(qū)域相吻合，證明了本文算法在尋找多個模體時仍具有較好的性能。

4 結(jié) 語

本文提出一種用于大規(guī)模ChIP?seq數(shù)據(jù)的替換顯露子串尋找問題的算法FastESE，通過測試集和控制集的比對找出出現(xiàn)頻率相對較高的子串，進一步搜索其（l，d）替換實例并進行聚類，篩選出具有較高信息量的子串即為替換顯露子串。使用老鼠胚胎干細(xì)胞的ChIP?seq數(shù)據(jù)進行實驗證明，F(xiàn)astESE可以在合理時間內(nèi)有效尋找到大規(guī)模ChIP?seq序列中的替換顯露子串。但是，在基因序列中這些替換顯露子串，即特殊的結(jié)合位點，往往具有不同功能，其功能性還需要進一步的實驗進行驗證。同時，在復(fù)雜的基因調(diào)控區(qū)域中具有特殊結(jié)構(gòu)與功能的顯露子串也是今后研究的重點之一。

參考文獻

[1] MARDIS E R. ChIP?seq： welcome to the new frontier [J]. Nat methods， 2007， 4（8）： 613?617.

[2] JIC C， CARSON MB， WANG Y， et al. A new exhaustive method and strategy for finding motifs in ChIP?enriched regions [J]. Plos one， 2014， 9（1）： e86044.

[3] JOHNSON D S. Genome?wide mapping of in vivo protein?DNA interactions [J]. Science， 2007， 316（5830）： 1497?1502.

[4] ZAMBELLI F， PESOLE G. Motif discovery and transcription factor binding sites before and after the next?generation sequencing era [J]. Briefings in bioinformatics， 2013， 14（2）： 225?237.

[5] PAVESI G， MAURI G. An algorithm for finding signals of unknown length in DNA sequences [J]. Bioinformatics， 2001， 17（Suppl 1）： S207?S214.

[6] HU M. On the detection and refinement of transcription factor binding sites using ChIP?Seq data [J]. Nucleic acids research， 2010， 38（7）： 2154?2167.

[7] REID J E， WERNISCH L. STEME： efficient em to find motifs in large data sets [J]. Nucleic acids research， 2011， 39（18）： 126?131.

[8] FISCHER J， HEUN V. Optimal string mining under frequency constraints， in knowledge discovery in databases [C]// Proceedings of 2006 PKDD. Heidelberg， Berlin： Springer， 2006：139?150.

[9] BULYK M L， JOHNSON P L. Nucleotides of transcription factor binding sites exert interdependent effects on the binding affinities of transcription factors [J]. Nucleic acids research， 2002， 30（5）： 1255?1261.

[10] SHAROV A A， Ko M S H. Exhaustive search for over?represented DNA sequence motifs with CisFinder [J]. DNA research， 2009， 16（5）： 261?273.

[11] VAN HEERINGEN S J， VEENSTRA G J C. Gimme motifs： a Denovo motif prediction pipeline for ChIP?sequencing experiments [J]. Bioinformatics， 2011， 27（2）： 270?271.

[12] YU Q. An efficient algorithm for discovering motifs in large DNA data sets [J]. IEEE transactions on nanobioscience， 2015， 14（5）： 535?44.

[13] THOMAS?CHOLLIER M. RSAT peak?motifs： motif analysis in full?size ChIP?seq datasets [J]. Nucleic acids research， 2012， 40（4）： e31?e31.

[14] SUN C， HUO H， YU Q， et al. An affinity propagation?based DNA motif discovery algorithm [J]. Biomed research international， 2015， 2015： 853461.