李 婷，張少強(qiáng)

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

隨著新技術(shù)的快速發(fā)展，基因組測(cè)序的成本下降，特別是轉(zhuǎn)錄因子的ChIP-seq技術(shù)的廣泛使用[1]，使得很多后生動(dòng)物和植物產(chǎn)生了海量的ChIP-seq數(shù)據(jù)集.盡管目前已有大量的預(yù)測(cè)順式調(diào)控元件和模塊的工具，但在大型基因組中，整合指數(shù)級(jí)增長的ChIP數(shù)據(jù)集，并在全基因組范圍預(yù)測(cè)順式調(diào)控元件和模塊，卻一直是具有挑戰(zhàn)性的計(jì)算問題[1-4].一定數(shù)量的轉(zhuǎn)錄因子常常組合起來，共同調(diào)控不同細(xì)胞類型、組織、發(fā)育階段和生理?xiàng)l件下的不同基因[5]，與這些共同調(diào)控的轉(zhuǎn)錄因子相結(jié)合的非編碼DNA位點(diǎn)(即順式調(diào)控元件)構(gòu)成了其順式調(diào)控模塊.大量的ChIP數(shù)據(jù)集中包含著一定的模塊組合信息，這些信息是由不同轉(zhuǎn)錄因子共同轉(zhuǎn)錄調(diào)控而形成的[6-7].因此，利用不同細(xì)胞類型、組織、發(fā)育階段和生理?xiàng)l件下的不同轉(zhuǎn)錄因子的大量ChIP數(shù)據(jù)集，就有可能通過對(duì)模體進(jìn)行整合以尋找共現(xiàn)模式，進(jìn)而對(duì)某種真核生物全基因組范圍順式調(diào)控模塊進(jìn)行從頭預(yù)測(cè).

本文基于果蠅已有的ChIP數(shù)據(jù)集，采用模體發(fā)現(xiàn)算法FisherNet及高性能并行的模體聚類算法CLIMP對(duì)果蠅的順式調(diào)控模塊進(jìn)行從頭預(yù)測(cè)，并與較新的DePCRM算法[8]進(jìn)行了比較.本文研究方法的流程如圖1所示.

圖1 順式調(diào)控模塊預(yù)測(cè)流程圖Fig.1 Flow chart of predicting CRMs

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

1.1 數(shù)據(jù)來源

由于果蠅常被用來研究動(dòng)物基因的轉(zhuǎn)錄調(diào)控，大量的順式調(diào)控元件和模塊已被實(shí)驗(yàn)驗(yàn)證，而且在過去的幾年中該生物已經(jīng)產(chǎn)生了大量的ChIP-chip和ChIP-seq數(shù)據(jù)，因此本文使用果蠅作為模式生物評(píng)估算法.為此，整理了來自56個(gè)不同轉(zhuǎn)錄因子的168個(gè)ChIP-chip和ChIP-seq數(shù)據(jù)集，這些數(shù)據(jù)集包含不同的發(fā)育階段(胚胎、幼蟲期1～3、蛹和雌雄成蟲)和不同實(shí)驗(yàn)條件下(熱休克等)的結(jié)果，其中：42個(gè)ChIP-chip和42 個(gè) ChIP-seq 數(shù)據(jù)集來自 modENCODE 項(xiàng)目[6，9]，38個(gè)ChIP-chip數(shù)據(jù)集來自Berkeley果蠅轉(zhuǎn)錄網(wǎng)絡(luò)項(xiàng)目(BDTNP)[10]，46個(gè)ChIP-chip數(shù)據(jù)集來自文獻(xiàn)[8].

1.2 數(shù)據(jù)預(yù)處理

利用peak-calling工具[11]查找ChIP數(shù)據(jù)中結(jié)合峰的序列，這些序列包含豐富的對(duì)應(yīng)轉(zhuǎn)錄因子的順式調(diào)控元件.將較短的結(jié)合峰從兩端延伸到3 000個(gè)堿基長的序列(這個(gè)長度與典型順式調(diào)控模塊的長度相當(dāng))，使得結(jié)合最高峰正好位于序列中部.除了ChIP實(shí)驗(yàn)的轉(zhuǎn)錄因子的順式調(diào)控元件外，擴(kuò)展的結(jié)合峰更可能包含輔助調(diào)控轉(zhuǎn)錄因子(在順式調(diào)控模塊中共同作用的轉(zhuǎn)錄因子)的順式調(diào)控元件.

2 算法步驟

數(shù)據(jù)預(yù)處理后的具體算法流程見圖2.

圖2 數(shù)據(jù)預(yù)處理后的具體算法流程圖Fig.2 Flow chart of detailed algorithm after data preprocessing

2.1 構(gòu)建模體相似多部圖

對(duì)于每組延伸后的結(jié)合峰序列數(shù)據(jù)集，運(yùn)用模體發(fā)現(xiàn)工具FisherNet算法[12]尋找大量的假定模體.對(duì)每個(gè)數(shù)據(jù)集輸出前k個(gè)最優(yōu)的模體，見圖2(a)，k默認(rèn)值為20.

對(duì)于預(yù)處理的每個(gè)數(shù)據(jù)集輸出的前20個(gè)最優(yōu)模體，以每個(gè)模體做為頂點(diǎn)，考慮到2個(gè)模體的頻率矩陣和位置權(quán)重矩陣，本文使用位置信息含量相似度量法 SPIC(similarity with position information contents)[13]計(jì)算不同數(shù)據(jù)集間模體的相似性(閾值為0.7)，SPIC度量法已被證實(shí)優(yōu)于其他度量公式[13]，若2個(gè)模體的相似度大于閾值，則連接2個(gè)模體，從而構(gòu)建模體相似多部圖，見圖2(b).數(shù)據(jù)集內(nèi)部模體之間不連邊，只計(jì)算不同數(shù)據(jù)集間模體的兩兩相似性.

構(gòu)建模體相似多部圖后，運(yùn)用雙向最佳匹配BDBM(bi-directional best match)算法尋找模體配對(duì)，見圖2(c)，其中，若一個(gè)模體與另外一個(gè)數(shù)據(jù)集中多個(gè)模體都最相似，則選取靠前的模體進(jìn)行配對(duì).

2.2 模體相似多部圖的CLIMP聚類

對(duì)于配對(duì)后的模體相似多部圖，運(yùn)用CLIMP算法[14]進(jìn)行團(tuán)(即每對(duì)頂點(diǎn)均連接的子圖)融合聚類，并形成聚類編號(hào)，見圖2(d).每個(gè)聚類中高度相似的模體分別來自于不同的數(shù)據(jù)集，這些相似的模體可能是同一轉(zhuǎn)錄因子在不同數(shù)據(jù)集的同一模體.因?yàn)橥晦D(zhuǎn)錄因子可能在多個(gè)ChIP數(shù)據(jù)集中作為輔調(diào)控因子或主調(diào)控因子出現(xiàn)，因此對(duì)應(yīng)的模體會(huì)在多個(gè)數(shù)據(jù)集中被反復(fù)識(shí)別.

2.3 構(gòu)建模體共現(xiàn)多部圖

對(duì)得到的團(tuán)融合聚類構(gòu)建模體共現(xiàn)多部圖，計(jì)算不同聚類中屬于相同數(shù)據(jù)集的每對(duì)模體的共現(xiàn)分?jǐn)?shù).對(duì)于數(shù)據(jù)集Md中的模體Md(i)和Md(j)，共現(xiàn)分?jǐn)?shù)Sc為

其中：|Md(i)|和|Md(j)|分別為模體Md(i)和Md(j)含有順式調(diào)控元件結(jié)合峰的數(shù)量；o(Md(i)，Md(j))代表這2個(gè)模體中都含有的順式調(diào)控元件的結(jié)合峰的數(shù)量.若共現(xiàn)分?jǐn)?shù)不小于閾值α，則視其為共現(xiàn)模體，將之連接，最終形成模體共現(xiàn)多部圖，見圖2(e).基于REDfly數(shù)據(jù)庫[15]已有順式調(diào)控模塊的訓(xùn)練，閾值α的取值為0.7.

2.4 模體共現(xiàn)多部圖的CLIMP聚類

對(duì)模體共現(xiàn)多部圖進(jìn)行CLIMP聚類，得到模塊類.聚類結(jié)果即為順式調(diào)控模塊，并按下式由小到大進(jìn)行排序

其中：M為聚類后的模塊；|M|為M中含有模體的數(shù)量；m為模塊中的模體；i(m)為模體m在團(tuán)融合聚類后的聚類編號(hào).SM的值越小，則順式調(diào)控模塊M就越可能是真實(shí)的.將少于2個(gè)模體的聚類舍棄.見圖2(f).

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

結(jié)合峰長度分布密度見圖3.圖中，虛線為結(jié)合峰長度分布密度，實(shí)線為結(jié)合峰長度的累積分布，可見結(jié)合峰的大部分長度約為1 000，有0.62%的結(jié)合峰長度大于5 000，由于其質(zhì)量不高，所以不使用這部分?jǐn)?shù)據(jù).由FisherNet查找的模體的信息含量分布密度見圖4.由圖4可見，162個(gè)數(shù)據(jù)集中的模體(有6個(gè)數(shù)據(jù)集包含模體少于2個(gè)，被丟棄)具有較高信息含量.在各個(gè)數(shù)據(jù)集輸出的前20個(gè)模體中，包含99個(gè)已知模體，并且被FisherNet程序優(yōu)先識(shí)別.

圖3 結(jié)合峰長度分布密度Fig.3 Distribution density of binding peak length

圖4 模體信息含量分布密度Fig.4 Distribution density of information content of motifs

將本算法(A)和DePCRM算法(B)應(yīng)用于162個(gè)ChIP數(shù)據(jù)集，模體和順式調(diào)控模塊預(yù)測(cè)結(jié)果見表1.其中，已知順式調(diào)控模塊數(shù)量為1 330個(gè)(REDfly數(shù)據(jù)庫).若一個(gè)已知的順式調(diào)控模塊與預(yù)測(cè)的順式調(diào)控模塊有至少一半長度是重疊的，則將其視為全覆蓋.

表1 本研究算法(A)和DePCRM算法(B)預(yù)測(cè)結(jié)果Tab.1 Predictions of algorithms of this research(A)and DePCRM(B)

由表1可見，在模體發(fā)現(xiàn)中，本算法輸出每個(gè)數(shù)據(jù)集中最優(yōu)的模體，得到了3 240個(gè)模體，其中包含1 214個(gè)已知的順式調(diào)控模塊(占已知數(shù)量的91.28%)；而DePCRM算法由于并未考慮模體的優(yōu)劣，因此輸出模體數(shù)量較多，為17890個(gè)，其中包含1 061個(gè)已知的順式調(diào)控模塊(占已知數(shù)量的79.77%).在順式調(diào)控模塊預(yù)測(cè)中，本算法得到的1 346個(gè)模塊中有1 103個(gè)已知模塊(占已知數(shù)量的82.93%)；而DePCRM算法得到的115 932個(gè)模塊中有947個(gè)已知模塊(占已知數(shù)量的71.20%).以上數(shù)據(jù)說明，本算法在順式調(diào)控模塊的預(yù)測(cè)中較DePCRM有更高的覆蓋率和敏感性.

順式調(diào)控模塊長度和相鄰順式調(diào)控元件間距離分布密度見圖 5(a)和(b).由圖 5(a)可見，本算法預(yù)測(cè)的順式調(diào)控模塊比已知的順式調(diào)控模塊的長度短.由圖5(b)可見，預(yù)測(cè)結(jié)果的相鄰順式調(diào)控元件間距離與已知的順式調(diào)控元件比較相似，一部分距離比已知的短.這表明可能遺漏了順式調(diào)控模塊中的某些順式調(diào)控元件，尤其是兩端的，這可能是由于ChIP數(shù)據(jù)沒有足夠多樣化的信息.

圖5 順式調(diào)控模塊長度預(yù)測(cè)結(jié)果Fig.5 Prediction results of CRM length

4 結(jié)論

本文利用大量的ChIP數(shù)據(jù)集實(shí)現(xiàn)了全基因組范圍的順式調(diào)控模塊的從頭預(yù)測(cè).通過識(shí)別最優(yōu)表達(dá)的、組合的模體，完成了對(duì)順式調(diào)控模塊的預(yù)測(cè).預(yù)測(cè)結(jié)果覆蓋了數(shù)據(jù)集中已知順式調(diào)控模塊的82.93%.這些預(yù)測(cè)的順式調(diào)控模塊比隨機(jī)選擇的序列更保守，更有可能具有調(diào)控功能.

與已有的DePCRM算法相比，本文采用了2個(gè)多部圖和2次CLIMP聚類，比DePCRM算法更簡(jiǎn)便快速.本算法不采用共現(xiàn)對(duì)的概念，克服了模體以偶數(shù)對(duì)出現(xiàn)的缺點(diǎn).當(dāng)有足夠多數(shù)量的、不同種類的其他真核生物ChIP數(shù)據(jù)集時(shí)，本算法可推廣到該類真核生物，用來預(yù)測(cè)其順式調(diào)控模塊.

參考文獻(xiàn)：

[1]PEPKE S，WOLD B，MORTAZAVI A.Computation for ChIP-seq and RNA-seq studies[J].Nature Methods，2009，6(11)：22-32.

[2]PARK P J.ChIP-seq：Advantages and challenges of a maturing technology[J].Nature Reviews Genetics，2009，10(10)：669-680.

[3]HAWKINS R D，HON G C，REN B.Next-generation genomics：An integrativeapproach[J].NatureReviewsGenetics，2010，11(7)：476-486.

[4]LAIRD P W.Principles and challenges of genome-wide DNA methylation analysis[J].Nature Reviews Genetics，2010，11(3)：191-203.

[5]MASTON G A，EVANS S K，GREEN M R.Transcriptional regulatory elements in the human genome[J].Annual Review of Genomics and Human Genetics，2006，7(1)：29-59.

[6]NEGRE N，BROWN C D，MA L J，et al.A cis-regulatory map of the drosophila genome[J].Nature，2011，471(7339)：527-531.

[7]GERSTEIN M B，LU Z J，NOSTRAND E L V，et al.Integrative analysis of the Caenorhabditis elegans genome by the modENCODE project[J].Science，2010，330(6012)：1775-1786.

[8]MENG N，TABARI E S，SU Z C.De novo prediction of cis-regulatory elements and modules through integrative analysis of a large number of ChIP datasets[J].BMC Genomics，2014，15(1)：1047-1066.

[9]CONSORTIUM T M，ROY S，ERNST J，et al.Identification of functional elements and regulatory circuits by Drosophila modENCODE[J].Science，2010，330(6012)：1787-1797.

[10]LI X Y，MACARTHUR S，BOURGON R，et al.Transcription factors bind thousands of active and inactive regions in the Drosophila blastoderm[J].Plos Biology，2008，6(2)：365-388.

[11]ZHANG Y，LIU T，MEYER C A，et al.Model-based analysis of ChIP-seq(MACS)[J].Genome Biology，2008，9(9)，DOI：10.1186/gb-2008-9-9-r137.

[12]張志紅.基于ChIP-seq數(shù)據(jù)集的順式調(diào)控模塊發(fā)現(xiàn)算法研究[D].天津：天津師范大學(xué)，2017.ZHANG Z H.Algorithm for Finding Cis-Regulatory Module Based on ChIP-seq Datasets[D].Tianjin：Tianjin Normal University，2017(in Chinese).

[13]ZHANG S Q，ZHOU X，et al.SPIC：A novel similarity metric for comparing transcription factor binding site motifs based on information contents[J].BMC Systems Biology，2013，7(2)：1-8.

[14]ZHANG S Q，CHEN Y.CLIMP：Clustering motifs via maximal cliques with parallel computing design[J].Plos One，2016，11(8)：1-17.

[15]IVAN A，HALFON M S，SINHA S.Computational discovery of cisregulatory modules in Drosophila，without prior knowledge of motifs[J].Genome Biology，2008，9(1)：1-17.