基于大量ChIP數據集的果蠅順式調控模塊的從頭預測

李 婷，張少強

(天津師范大學 計算機與信息工程學院，天津 300387)

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

本文基于果蠅已有的ChIP數據集，采用模體發現算法FisherNet及高性能并行的模體聚類算法CLIMP對果蠅的順式調控模塊進行從頭預測，并與較新的DePCRM算法[8]進行了比較.本文研究方法的流程如圖1所示.

圖1 順式調控模塊預測流程圖Fig.1 Flow chart of predicting CRMs

1 數據來源與預處理

1.1 數據來源

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

1.2 數據預處理

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

2 算法步驟

數據預處理后的具體算法流程見圖2.

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

2.1 構建模體相似多部圖

對于每組延伸后的結合峰序列數據集，運用模體發現工具FisherNet算法[12]尋找大量的假定模體.對每個數據集輸出前k個最優的模體，見圖2(a)，k默認值為20.

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

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

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

對于配對后的模體相似多部圖，運用CLIMP算法[14]進行團(即每對頂點均連接的子圖)融合聚類，并形成聚類編號，見圖2(d).每個聚類中高度相似的模體分別來自于不同的數據集，這些相似的模體可能是同一轉錄因子在不同數據集的同一模體.因為同一轉錄因子可能在多個ChIP數據集中作為輔調控因子或主調控因子出現，因此對應的模體會在多個數據集中被反復識別.

2.3 構建模體共現多部圖

對得到的團融合聚類構建模體共現多部圖，計算不同聚類中屬于相同數據集的每對模體的共現分數.對于數據集Md中的模體Md(i)和Md(j)，共現分數Sc為

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

2.4 模體共現多部圖的CLIMP聚類

對模體共現多部圖進行CLIMP聚類，得到模塊類.聚類結果即為順式調控模塊，并按下式由小到大進行排序

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

3 實驗結果

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

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

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

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

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

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

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

圖5 順式調控模塊長度預測結果Fig.5 Prediction results of CRM length

4 結論

本文利用大量的ChIP數據集實現了全基因組范圍的順式調控模塊的從頭預測.通過識別最優表達的、組合的模體，完成了對順式調控模塊的預測.預測結果覆蓋了數據集中已知順式調控模塊的82.93%.這些預測的順式調控模塊比隨機選擇的序列更保守，更有可能具有調控功能.

與已有的DePCRM算法相比，本文采用了2個多部圖和2次CLIMP聚類，比DePCRM算法更簡便快速.本算法不采用共現對的概念，克服了模體以偶數對出現的缺點.當有足夠多數量的、不同種類的其他真核生物ChIP數據集時，本算法可推廣到該類真核生物，用來預測其順式調控模塊.

參考文獻：

[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數據集的順式調控模塊發現算法研究[D].天津：天津師范大學，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.