基于絕緣密度HiCID算法的染色質(zhì)接觸域邊界檢測

黃月月，豐繼華，劉 珂，范力棟

（云南民族大學(xué)電氣信息工程學(xué)院，昆明650504）

基因組的空間結(jié)構(gòu)與其生物功能密切相關(guān)。高通量測序Hi?C 技術(shù)的出現(xiàn)使得研究全基因組成對基因座之間的交互接觸作用成為可能，為完整勾畫出染色體在三維空間上的交互網(wǎng)絡(luò)奠定了基礎(chǔ)［1］。對酵母、果蠅、小鼠和人類染色體的研究發(fā)現(xiàn)［2?6］，染色體是由數(shù)百kb 到數(shù)十Mb 大小不等的染色質(zhì)接觸域構(gòu)成，從而形成染色質(zhì)區(qū)室、拓?fù)潢P(guān)聯(lián)域（Topologically associat?ing domain，TAD）和染色體環(huán)等基本結(jié)構(gòu)。高分辨率Hi?C 交互作用圖譜揭示了接觸域內(nèi)部的相互作用強(qiáng)，不同接觸域間的相互作用弱［7］。

高通量生物學(xué)數(shù)據(jù)龐大且復(fù)雜，交互測量過程中產(chǎn)生的噪聲會對不同類型的生物網(wǎng)絡(luò)造成干擾，如Hi?C 交互網(wǎng)絡(luò)、細(xì)胞間交互網(wǎng)絡(luò)［8］和PPI 網(wǎng)絡(luò)［9］等。噪聲的存在尤其對下游性能的分析產(chǎn)生不利的影響，因此首先需借助基于網(wǎng)絡(luò)的降噪方法對原始交互數(shù)據(jù)進(jìn)行預(yù)處理。網(wǎng)絡(luò)增強(qiáng)NE［10］可有效改善生物網(wǎng)絡(luò)的質(zhì)量，目前已將NE 應(yīng)用于人類組織的22 個基因相互作用網(wǎng)絡(luò)的基因功能預(yù)測、改善Hi?C 網(wǎng)絡(luò)以促進(jìn)拓?fù)溆虻淖R別以及提高細(xì)粒度物種鑒定的準(zhǔn)確度等研究領(lǐng)域。全基因組研究發(fā)現(xiàn)，染色質(zhì)接觸域邊界除了富集CTCF的結(jié)合位點和組蛋白化學(xué)修飾外，還有大量的管家基因、tRNAs、SINE 反轉(zhuǎn)錄轉(zhuǎn)座子等DNA 元件［11］。因此利用Hi?C 數(shù)據(jù)進(jìn)行接觸域邊界檢測與定位，對基因調(diào)控、基因組相互作用以及基因功能等方面的研究具有十分重要的生物學(xué)意義。目前，染色質(zhì)接觸域的檢測識別已涌現(xiàn)出許多檢測算法。根據(jù)檢測原理的不同，主要將接觸域檢測模型分為兩類：基于一維統(tǒng)計量的檢測模型和基于二維接觸矩陣的檢測模型?；谝痪S統(tǒng)計量的檢測模型：Dixon 等［12］引入方向性指數(shù)（Di?rectionality index，DI）和隱馬爾可夫模型（Hidden mar?kov model，HMM）來推斷基因組中接觸域的位置；Shin 等［13］提出TopDom 算法檢測接觸域及其邊界；Chen 等［14］在TopDom 算法基礎(chǔ)上提出HiCDB 算法，進(jìn)而識別多分辨率的Hi?C 數(shù)據(jù)的接觸域邊界。基于二維接觸矩陣的算法種類比較繁多，檢測原理各異，目前較為常用的算法有ClusterTAD、HiCSeg、IC?Finder 和Arrowhead等［15?18］。本文在現(xiàn)有的HiCDB算法基礎(chǔ)上，提出邊界識別精度更高且性能更好的絕緣密度算法HiCID。

1 材料與方法

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

實驗采用分辨率40 kb 的小鼠mESC（胚胎干細(xì)胞）、mCO（皮質(zhì)細(xì)胞）細(xì)胞系Hi?C 數(shù)據(jù)（http：//chromo?some.sdsc.edu/mouse/hi?c/download.html）和分辨率為10 kb 的人類IMR90（胚肺成纖維）細(xì)胞系Hi?C 數(shù)據(jù)（NCBI 登錄號：GSE35156）。不同細(xì)胞系對應(yīng)的CTCF和各種組蛋白修飾ChIP?seq數(shù)據(jù)均可以從ENCODE 數(shù)據(jù)庫下載。此外，分別選擇人類hg19 和小鼠細(xì)胞系mm10 作為參考基因組。網(wǎng)絡(luò)增強(qiáng)實驗代碼鏈接（http：//snap.stanford.edu/ne/）；HiCDB 程 序 代 碼（https：//github.com/ChenFengling/HiCDB）。

1.2 HiCID算法

1.2.1 網(wǎng)絡(luò)增強(qiáng)去噪

網(wǎng)絡(luò)增強(qiáng)（Network enhancement，NE）算法是一種基于動態(tài)擴(kuò)散過程的網(wǎng)絡(luò)去噪算法，輸入是被噪聲干擾的無向加權(quán)網(wǎng)絡(luò)，通過擴(kuò)散過程不斷更新網(wǎng)絡(luò)連接權(quán)重，直至權(quán)重收斂生成一個與輸入網(wǎng)絡(luò)具有相同節(jié)點數(shù)的新網(wǎng)絡(luò)結(jié)構(gòu)。NE 算法的核心是利用對稱、半正定、且能誘導(dǎo)稀疏性的雙隨機(jī)矩陣（Doubly stochastic matrix，DSM）算子來增強(qiáng)網(wǎng)絡(luò)節(jié)點之間的相似性，進(jìn)而提高網(wǎng)絡(luò)信噪比、改善下游網(wǎng)絡(luò)分析的性能。

1.2.2 絕緣密度（Insulation density）

定義的統(tǒng)計量絕緣密度（Insulation density，ID）可以簡單有效地計算每個基因位點在w 窗長范圍內(nèi)絕緣強(qiáng)度的密度分布，它將二維的Hi?C 數(shù)據(jù)矩陣轉(zhuǎn)換為一維矢量，為了克服偶然因素對實驗結(jié)果造成的影響，使識別的接觸域邊界更具說服力，在不同大小的窗口下計算平均絕緣密度（Average insulation density，AID）：

其中，m表示不同大小窗口的數(shù)量，第m個窗口長度為wm，s（k，k+1）為基因位點的中心，U、D 和B 分別為s上、下游和之間區(qū)域窗口w內(nèi)所有Hi?C交互頻率之和。

最佳窗口數(shù)量m 的確定可以借助接觸域質(zhì)量TADquality 來約束，在最佳窗口數(shù)量下，接觸域獲得最高質(zhì)量分?jǐn)?shù)，此時接觸域內(nèi)各交互頻率的相似性達(dá)到最大化。拓?fù)溆蛸|(zhì)量公式如下：

其中intra（i）表示染色體上第i個接觸域內(nèi)的平均接觸頻率，intra（i，j）表示位于第i個接觸域和第j個接觸域之間區(qū)域的接觸頻率平均值。

1.2.3 去除背景噪聲及局部極值檢測

AID 信號需要進(jìn)一步去除背景噪聲以增強(qiáng)信號穩(wěn)定性。通過線性擬合AID 信號的局部最小值生成雙層下包絡(luò)，AID信號減去雙層下包絡(luò)，最終獲得平滑背景的局部絕緣密度（Local insulation density，LID）。AID信號的局部極大峰值可通過MATLAB 內(nèi)置函數(shù)find?peaks 或自動多尺度峰值檢測［19］（Automatic multiscale?based peak detection，AMPD）算法進(jìn)行識別，峰值位置對應(yīng)的基因位點即為候選接觸域邊界。

1.2.4 局部極值的截止閾值

利用CTCF［20］與其他組蛋白標(biāo)記的富集信息共同確定閾值：首先對LID 的局部極大峰值按降序進(jìn)行排序，如果LID 峰值所在si，j位置處有CTCF、H3K4me1 和H3K4me3等出現(xiàn)，則注釋為1，否則為0，由此獲得各分向量S：

如果H 表示富集CTCF、H3K4me1 或H3K4me3 等元素的LID峰值集合，則：

其中，i ∈{1 ,2,…,h}，j ∈{1 ,2,…,n}，h表示CTCF、H3K4me1或H3K4me3 等元素種類的個數(shù)，n 表示所有LID 峰值的個數(shù)，nhit表示有CTCF、H3K4me1、H3K4me3 等存在的LID峰值個數(shù)。

接著利用概率知識來計算豐度（Enrichment score，ES），平均豐度取得最大值時的LID 設(shè)為局部極大峰值的截止閾值，而被濾除的局部極大峰值可能是實驗噪聲等不確定因素引起的，這一定程度上有助于提升邊界識別的準(zhǔn)確性。

H 中第i 個局部極大峰值LIDi出現(xiàn)的概率為；由于n 為識別的接觸域邊界的總數(shù)，因此表示沒有CTCF 或其他組蛋白標(biāo)記存在時每個峰值出現(xiàn)的概率。

1.2.5 方法流程

接觸域邊界檢測方法流程見圖1。

2 結(jié)果與分析

2.1 網(wǎng)絡(luò)增強(qiáng)去噪效果對比

圖2 是在人類H1 細(xì)胞系chr1 染色體上截取的［100，300］基因區(qū)間內(nèi)的Hi?C 圖譜。從圖2 可以直觀地看出，原始Hi?C 數(shù)據(jù)生成的熱圖被噪聲嚴(yán)重干擾，網(wǎng)絡(luò)節(jié)點間的連接強(qiáng)度較弱，接觸域分層結(jié)構(gòu)模糊，相鄰接觸域之間的相互連接邊界不清晰。經(jīng)降噪處理后的Hi?C 熱圖對角線上的連續(xù)方形區(qū)域顏色更突出，而且不同區(qū)域之間的邊界輪廓更清晰，域內(nèi)嵌套的子TAD結(jié)構(gòu)更容易識別。

圖1 接觸域邊界檢測流程圖Figure 1 The flowchart of topology domain boundary detection

圖2 網(wǎng)絡(luò)增強(qiáng)前后效果對比圖Figure 2 The comparison between before and after network enhancement

2.2 不同的絕緣系數(shù)和峰值檢測效果對比

為了進(jìn)一步比較不同絕緣系數(shù)或峰值檢測方法之間的差異性，在圖1 基因區(qū)間內(nèi)分別繪制由相對絕緣系數(shù)RI 獲得的LRI 信號和由絕緣密度公式ID 獲得的LID 信號。經(jīng)對比發(fā)現(xiàn)，LRI 與LID 的曲線走勢基本一致，但由絕緣密度公式獲得的LID 絕緣性更高，峰值銳化更明顯，同時能識別被相對絕緣公式遺漏的少量峰值。AMPD 識別的峰值與findpeaks 函數(shù)的檢測結(jié)果基本一致（圖3）。

圖3 不同的絕緣系數(shù)和峰值檢測效果對比Figure 3 The comparison of different insulation coefficients and peak detection methods

2.3 TADquality確定窗口大小和數(shù)量

針對人類H1 細(xì)胞系的22 條常染色體，分別在不同窗口長度下計算所有接觸域內(nèi)的平均TADquality。從圖4 可以看出，隨著窗口長度的不斷增大，TADqual?ity 在數(shù)值上隨之增加。當(dāng)窗口達(dá)到6 個bin 的長度以后，TADquality 逐漸有收斂的趨勢；當(dāng)w=7 時，獲得了最高的TADquality。多個窗口下的平均TADquality 與單個窗口的TADquality 變化趨勢基本一致，多個窗口取平均考慮了Hi?C 數(shù)據(jù)在不同窗口下的伸縮性，可以有效地減少實驗誤差。

圖4 TADquality與窗口長度的關(guān)系圖Figure 4 The relationship between TADquality and window length

2.4 GSEA確定截止閾值

圖5 展示了H1 細(xì)胞 系chr1 染色體上CTCF 和8 種常見的組蛋白修飾在接觸域邊界附近的富集狀態(tài)。以CTCF 為例，將CTCF 在基因組范圍內(nèi)的峰值分布信息作為輸入，以單個域邊界為中心，分別統(tǒng)計上、下游50 個bin 區(qū)間內(nèi)CTCF 峰值出現(xiàn)的頻率。由于域邊界數(shù)量較多，則需要計算每個bin 內(nèi)平均峰值個數(shù)來表征CTCF 的平均富集情況。從圖5 可以看出：邊界處出現(xiàn)CTCF 峰值的頻率較高，說明CTCF 在域邊界處有較強(qiáng)程度的富集；H3K4me2、H3K4me3、H4K20me1、H3K36me3 以及H3K9ac 等組蛋白在邊界處的平均峰值數(shù)量顯著增多，它們通常分布在轉(zhuǎn)錄起始位點附近的啟動子區(qū)域，用于激活基因的表達(dá)；H3K27ac 和H3K4me1 共同作為活性基因增強(qiáng)子的標(biāo)志，在邊界附近的基因區(qū)間內(nèi)有較為平穩(wěn)的富集狀態(tài)；而H3K27me3 和H3K9me3 與基因抑制有關(guān)，它們在邊界處的峰值變化不明顯，甚至有低谷（或損耗）出現(xiàn)。

圖5 接觸域邊界附近CTCF和8種組蛋白峰值富集狀態(tài)Figure 5 The peak enrichment status of CTCF and proteins near the topological domain boundary

盡管大多數(shù)邊界富集了CTCF結(jié)合位點，但單一變量CTCF 還不足以識別域邊界。以基因組峰值分布信息為例，人類H1 細(xì)胞系中的22 條常染色體上共有CTCF 峰值63 863 個，只有大約30.26%的CTCF 峰值位點位于接觸域邊界的附近區(qū)域，這種相對較低的數(shù)據(jù)利用率可能與識別接觸域時設(shè)置的嚴(yán)格閾值有關(guān)。然而，CTCF 不存在的基因位點處也可能會有接觸域存在，因此需要借助其他協(xié)變量來完成截止閾值的設(shè)定。在邊界處有顯著變化（明顯增加或損耗）的組蛋白共有6 種，其中CTCF 和H4K20me1 的組合最高，將CTCF的利用率提高到了35.25%，而CTCF和H4K9me3也在一定程度上將利用率提高到了33.05%，說明缺乏H3K9me3也可以作為識別域邊界的重要特征（圖6）。

圖6 CTCF與各組蛋白修飾組合對數(shù)據(jù)利用率的影響Figure 6 The effect of CTCF and histone modification combinations on data utilization

在40 kb 的人類H1 細(xì)胞系中，分別采用不同的協(xié)變量組合方式確定截止閾值，進(jìn)而調(diào)用作用域的邊界。表1顯示，CTCF作為唯一變量識別的接觸域邊界數(shù)量較多，但只有65.28%的邊界是保守性的。一般地，保守性越高說明域邊界的預(yù)測結(jié)果越準(zhǔn)確。CTCF 與其他組蛋白的結(jié)合會提高截止閾值，使得域邊界數(shù)量減少，但保守性邊界比重增加，表明借助組蛋白修飾的確可以有效提升預(yù)測結(jié)果的準(zhǔn)確性。另外，由于染色體的長度固定，所以相鄰域邊界之間的距離會隨著邊界數(shù)量的減少而增加。通常兩個相鄰邊界之間的距離越小，說明域結(jié)構(gòu)越精細(xì)，能捕捉到嵌套TAD 的概率就越大。

表1 組蛋白修飾對域邊界檢測結(jié)果的影響Table 1 The effect of histone modifications on domain boundary detection

利用類似GSEA 的算法確定截止閾值。圖7（a）是經(jīng)峰值檢測識別的所有局部極大峰值的降序排列，圖7（b）為降序排列后的峰值位點的富集分?jǐn)?shù)，不同的協(xié)變量組合方式對應(yīng)的富集分?jǐn)?shù)的極大值不同。

圖7 GSEA算法確定截止閾值Figure 7 The determination of cut?off threshold by GSEA

2.5 性能評價指標(biāo)

2.5.1 不同域邊界檢測方法的一致性比較

圖8 是HiCID 算法與其他3 種接觸域檢測算法（HiCDB、TopDom 和DomainCaller［22］）的性能比對分析。表2 和表3 分別統(tǒng)計了人類和小鼠不同細(xì)胞系采用不同域邊界檢測算法的一致性檢測結(jié)果，以40 kb人類IMR90 為例，由于CTCF 與組蛋白H3K9me3 的組合最高將保守性提高到87.64%，所以該組合被應(yīng)用于HiCID 算法中以篩選域邊界。與hESC 細(xì)胞系中的檢測結(jié)果相似，相較于原HiCDB，NE 降噪后的HiCDB 算法（HiCDB+NE）識別的邊界數(shù)量有所減少，但邊界保守性基本保持穩(wěn)定。HiCID 算法識別的域邊界數(shù)量介于HiCDB 算法和HiCDB+NE 之間，但保守性邊界占邊界總數(shù)的百分比均高于此兩種算法。盡管TopDom 和DomainCaller識別的域邊界數(shù)量有限，但保守性邊界居多，而且針對不同的物種有較為穩(wěn)定的檢測性能，但是這兩種算法計算的相鄰邊界間的平均距離較大，不能識別更加精細(xì)的接觸域結(jié)構(gòu)。同時發(fā)現(xiàn)，HiCID 算法應(yīng)用于10 kb 的IMR90 細(xì)胞系時，與其他3 種算法相比，不僅識別的邊界數(shù)量最多、保守性最高，而且相鄰邊界之間的距離最接近。由此推測，HiCID 算法對較高分辨率的Hi?C數(shù)據(jù)更有效。

圖8 人類hESC和IMR90細(xì)胞系中不同域邊界檢測方法的一致性比較Figure 8 The comparison of consistency methods for different domains of human hESC and IMR90

表2 人類細(xì)胞系不同域邊界檢測方法的一致性比較Table 2 The comparison of consistency methods for different domains of human hESC and IMR90

表3 小鼠細(xì)胞系不同域邊界檢測算法一致性比較Table 3 Comparison of the consistency of different cell boundary detection algorithms in mouse cell lines

2.5.2 域邊界檢測結(jié)果的準(zhǔn)確性

高活性的組蛋白修飾可用來標(biāo)記作用域的起點和終止位點，即域邊界位置。因此，通過統(tǒng)計比較邊界附近各組蛋白修飾的數(shù)量可進(jìn)一步判斷域邊界的準(zhǔn)確性。圖9 顯示，在整個統(tǒng)計區(qū)間內(nèi)，對于4 種不同的組蛋白修飾，表征HiCID 算法的藍(lán)色曲線與表征HiCDB+NE 算法的紅色曲線在個別區(qū)間內(nèi)有交錯重疊，但藍(lán)色曲線的整體幅值略高于紅色曲線的幅值。

為了從數(shù)值上體現(xiàn)兩種算法在性能上的差異，分別在hESC 細(xì)胞系的22 條常染色體和小鼠Cortex 細(xì)胞系19 條常染色體上計算了邊界附近單位bin 內(nèi)平均峰值的個數(shù)。從圖10 可以看出，HiCID 算法檢測的5 種組蛋白的平均峰值的個數(shù)均大于HiCDB+NE 算法檢測的平均峰值個數(shù)，說明在識別域邊界時，HiCID 算法在整個基因組都普遍有效。

圖9 4種組蛋白修飾在hESC域邊界附近平均峰值個數(shù)Figure 9 Average number of peaks of four histone modifications near the boundary of hESC domain

2.5.3 接觸域邊界檢測結(jié)果

圖11是HiCDB與HiCID識別拓?fù)溆蜻吔绲慕Y(jié)果對比圖。人類hESC 細(xì)胞系chr18 染色體的［67 100 000，71 100 000］基因區(qū)間約有154個bin（分辨率為40 kb），Hi?C熱圖中用淺藍(lán)色的點標(biāo)注了兩種算法識別的域邊界位置，而深藍(lán)色的點表示在21 種人類細(xì)胞系中都存在的保守域邊界。從圖11 可以直觀地看出，改進(jìn)前后HiCDB 識別的域邊界有部分重合，尤其對保守性邊界標(biāo)注基本一致，而HiCID不僅能識別邊界輪廓清晰的域結(jié)構(gòu)，對域內(nèi)出現(xiàn)的嵌套TAD也能有效識別。

3 討論與結(jié)論

當(dāng)前表觀遺傳學(xué)領(lǐng)域涌現(xiàn)出眾多域檢測算法，這些方法為檢測拓?fù)溆蚣捌溥吔缣峁┝烁嗟倪x擇，但大多數(shù)檢測方法在應(yīng)用過程中都有一定的局限性，比如：有的算法需要對多個參數(shù)進(jìn)行調(diào)整；有的不適用于有多種分辨率的Hi?C數(shù)據(jù)；或者算法的可重復(fù)性差、運行時間成本較高、檢測結(jié)果的準(zhǔn)確性有待提高。

圖10 hESC和Cortex細(xì)胞系各組蛋白平均出現(xiàn)的峰值個數(shù)Figure 10 The average number of peaks of each histone on 19 chromosomes in mice Cortex and hESC

圖11 HiCDB與HiCID接觸域邊界檢測結(jié)果Figure 11 The detection results of HiCDB and HiCID contact domain boundaries

與代表性的HiCDB 算法相比較，本文提出的Hi?CID算法利用網(wǎng)絡(luò)增強(qiáng)技術(shù)對原始Hi?C數(shù)據(jù)進(jìn)行降噪預(yù)處理，構(gòu)建絕緣密度統(tǒng)計量確定邊界特征，具有不易漏檢和冗余性好的特點。針對HiCDB算法在篩選閾值時只考慮單一協(xié)變量CTCF，本文在識別過程中增加了組蛋白修飾信息，進(jìn)一步提升了算法穩(wěn)健性。此外，HiCID 適用范圍廣泛，對于不同分辨率下的Hi?C數(shù)據(jù)，最優(yōu)窗口參數(shù)的選擇可以利用TADquality有效約束。

由于現(xiàn)有實驗數(shù)據(jù)的局限性，本文提出的算法對邊界元件信息的融合程度還不夠，如管家基因、tRNAs、SINE 反轉(zhuǎn)錄轉(zhuǎn)座子等DNA 元件的分布信息未能引入算法，因此在處理不同物種或不同細(xì)胞系時算法性能還不夠穩(wěn)定。此外，除算法本身的局限性之外，目前關(guān)于染色質(zhì)相互作用域的大小和劃分還沒有一個標(biāo)準(zhǔn)的定義，在同一個作用域中會包含大量嵌套的TAD。此外，在細(xì)胞過程中會發(fā)生改變和重組，給有效定位域邊界帶來了挑戰(zhàn)。