基于绝缘密度HiCID算法的染色质接触域边界检测

黄月月，丰继华，刘 珂，范力栋

（云南民族大学电气信息工程学院，昆明650504）

基因组的空间结构与其生物功能密切相关。高通量测序Hi⁃C 技术的出现使得研究全基因组成对基因座之间的交互接触作用成为可能，为完整勾画出染色体在三维空间上的交互网络奠定了基础［1］。对酵母、果蝇、小鼠和人类染色体的研究发现［2⁃6］，染色体是由数百kb 到数十Mb 大小不等的染色质接触域构成，从而形成染色质区室、拓扑关联域（Topologically associat⁃ing domain，TAD）和染色体环等基本结构。高分辨率Hi⁃C 交互作用图谱揭示了接触域内部的相互作用强，不同接触域间的相互作用弱［7］。

高通量生物学数据庞大且复杂，交互测量过程中产生的噪声会对不同类型的生物网络造成干扰，如Hi⁃C 交互网络、细胞间交互网络［8］和PPI 网络［9］等。噪声的存在尤其对下游性能的分析产生不利的影响，因此首先需借助基于网络的降噪方法对原始交互数据进行预处理。网络增强NE［10］可有效改善生物网络的质量，目前已将NE 应用于人类组织的22 个基因相互作用网络的基因功能预测、改善Hi⁃C 网络以促进拓扑域的识别以及提高细粒度物种鉴定的准确度等研究领域。全基因组研究发现，染色质接触域边界除了富集CTCF的结合位点和组蛋白化学修饰外，还有大量的管家基因、tRNAs、SINE 反转录转座子等DNA 元件［11］。因此利用Hi⁃C 数据进行接触域边界检测与定位，对基因调控、基因组相互作用以及基因功能等方面的研究具有十分重要的生物学意义。目前，染色质接触域的检测识别已涌现出许多检测算法。根据检测原理的不同，主要将接触域检测模型分为两类：基于一维统计量的检测模型和基于二维接触矩阵的检测模型。基于一维统计量的检测模型：Dixon 等［12］引入方向性指数（Di⁃rectionality index，DI）和隐马尔可夫模型（Hidden mar⁃kov model，HMM）来推断基因组中接触域的位置；Shin 等［13］提出TopDom 算法检测接触域及其边界；Chen 等［14］在TopDom 算法基础上提出HiCDB 算法，进而识别多分辨率的Hi⁃C 数据的接触域边界。基于二维接触矩阵的算法种类比较繁多，检测原理各异，目前较为常用的算法有ClusterTAD、HiCSeg、IC⁃Finder 和Arrowhead等［15⁃18］。本文在现有的HiCDB算法基础上，提出边界识别精度更高且性能更好的绝缘密度算法HiCID。

1 材料与方法

1.1 数据来源

实验采用分辨率40 kb 的小鼠mESC（胚胎干细胞）、mCO（皮质细胞）细胞系Hi⁃C 数据（http：//chromo⁃some.sdsc.edu/mouse/hi⁃c/download.html）和分辨率为10 kb 的人类IMR90（胚肺成纤维）细胞系Hi⁃C 数据（NCBI 登录号：GSE35156）。不同细胞系对应的CTCF和各种组蛋白修饰ChIP⁃seq数据均可以从ENCODE 数据库下载。此外，分别选择人类hg19 和小鼠细胞系mm10 作为参考基因组。网络增强实验代码链接（http：//snap.stanford.edu/ne/）；HiCDB 程 序 代 码（https：//github.com/ChenFengling/HiCDB）。

1.2 HiCID算法

1.2.1 网络增强去噪

网络增强（Network enhancement，NE）算法是一种基于动态扩散过程的网络去噪算法，输入是被噪声干扰的无向加权网络，通过扩散过程不断更新网络连接权重，直至权重收敛生成一个与输入网络具有相同节点数的新网络结构。NE 算法的核心是利用对称、半正定、且能诱导稀疏性的双随机矩阵（Doubly stochastic matrix，DSM）算子来增强网络节点之间的相似性，进而提高网络信噪比、改善下游网络分析的性能。

1.2.2 绝缘密度（Insulation density）

定义的统计量绝缘密度（Insulation density，ID）可以简单有效地计算每个基因位点在w 窗长范围内绝缘强度的密度分布，它将二维的Hi⁃C 数据矩阵转换为一维矢量，为了克服偶然因素对实验结果造成的影响，使识别的接触域边界更具说服力，在不同大小的窗口下计算平均绝缘密度（Average insulation density，AID）：

其中，m表示不同大小窗口的数量，第m个窗口长度为wm，s（k，k+1）为基因位点的中心，U、D 和B 分别为s上、下游和之间区域窗口w内所有Hi⁃C交互频率之和。

最佳窗口数量m 的确定可以借助接触域质量TADquality 来约束，在最佳窗口数量下，接触域获得最高质量分数，此时接触域内各交互频率的相似性达到最大化。拓扑域质量公式如下：

其中intra（i）表示染色体上第i个接触域内的平均接触频率，intra（i，j）表示位于第i个接触域和第j个接触域之间区域的接触频率平均值。

1.2.3 去除背景噪声及局部极值检测

AID 信号需要进一步去除背景噪声以增强信号稳定性。通过线性拟合AID 信号的局部最小值生成双层下包络，AID信号减去双层下包络，最终获得平滑背景的局部绝缘密度（Local insulation density，LID）。AID信号的局部极大峰值可通过MATLAB 内置函数find⁃peaks 或自动多尺度峰值检测［19］（Automatic multiscale⁃based peak detection，AMPD）算法进行识别，峰值位置对应的基因位点即为候选接触域边界。

1.2.4 局部极值的截止阈值

利用CTCF［20］与其他组蛋白标记的富集信息共同确定阈值：首先对LID 的局部极大峰值按降序进行排序，如果LID 峰值所在si，j位置处有CTCF、H3K4me1 和H3K4me3等出现，则注释为1，否则为0，由此获得各分向量S：

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

其中，i ∈{1 ,2,…,h}，j ∈{1 ,2,…,n}，h表示CTCF、H3K4me1或H3K4me3 等元素种类的个数，n 表示所有LID 峰值的个数，nhit表示有CTCF、H3K4me1、H3K4me3 等存在的LID峰值个数。

接着利用概率知识来计算丰度（Enrichment score，ES），平均丰度取得最大值时的LID 设为局部极大峰值的截止阈值，而被滤除的局部极大峰值可能是实验噪声等不确定因素引起的，这一定程度上有助于提升边界识别的准确性。

H 中第i 个局部极大峰值LIDi出现的概率为；由于n 为识别的接触域边界的总数，因此表示没有CTCF 或其他组蛋白标记存在时每个峰值出现的概率。

1.2.5 方法流程

接触域边界检测方法流程见图1。

2 结果与分析

2.1 网络增强去噪效果对比

图2 是在人类H1 细胞系chr1 染色体上截取的［100，300］基因区间内的Hi⁃C 图谱。从图2 可以直观地看出，原始Hi⁃C 数据生成的热图被噪声严重干扰，网络节点间的连接强度较弱，接触域分层结构模糊，相邻接触域之间的相互连接边界不清晰。经降噪处理后的Hi⁃C 热图对角线上的连续方形区域颜色更突出，而且不同区域之间的边界轮廓更清晰，域内嵌套的子TAD结构更容易识别。

图1 接触域边界检测流程图Figure 1 The flowchart of topology domain boundary detection

图2 网络增强前后效果对比图Figure 2 The comparison between before and after network enhancement

2.2 不同的绝缘系数和峰值检测效果对比

为了进一步比较不同绝缘系数或峰值检测方法之间的差异性，在图1 基因区间内分别绘制由相对绝缘系数RI 获得的LRI 信号和由绝缘密度公式ID 获得的LID 信号。经对比发现，LRI 与LID 的曲线走势基本一致，但由绝缘密度公式获得的LID 绝缘性更高，峰值锐化更明显，同时能识别被相对绝缘公式遗漏的少量峰值。AMPD 识别的峰值与findpeaks 函数的检测结果基本一致（图3）。

图3 不同的绝缘系数和峰值检测效果对比Figure 3 The comparison of different insulation coefficients and peak detection methods

2.3 TADquality确定窗口大小和数量

针对人类H1 细胞系的22 条常染色体，分别在不同窗口长度下计算所有接触域内的平均TADquality。从图4 可以看出，随着窗口长度的不断增大，TADqual⁃ity 在数值上随之增加。当窗口达到6 个bin 的长度以后，TADquality 逐渐有收敛的趋势；当w=7 时，获得了最高的TADquality。多个窗口下的平均TADquality 与单个窗口的TADquality 变化趋势基本一致，多个窗口取平均考虑了Hi⁃C 数据在不同窗口下的伸缩性，可以有效地减少实验误差。

图4 TADquality与窗口长度的关系图Figure 4 The relationship between TADquality and window length

2.4 GSEA确定截止阈值

图5 展示了H1 细胞 系chr1 染色体上CTCF 和8 种常见的组蛋白修饰在接触域边界附近的富集状态。以CTCF 为例，将CTCF 在基因组范围内的峰值分布信息作为输入，以单个域边界为中心，分别统计上、下游50 个bin 区间内CTCF 峰值出现的频率。由于域边界数量较多，则需要计算每个bin 内平均峰值个数来表征CTCF 的平均富集情况。从图5 可以看出：边界处出现CTCF 峰值的频率较高，说明CTCF 在域边界处有较强程度的富集；H3K4me2、H3K4me3、H4K20me1、H3K36me3 以及H3K9ac 等组蛋白在边界处的平均峰值数量显著增多，它们通常分布在转录起始位点附近的启动子区域，用于激活基因的表达；H3K27ac 和H3K4me1 共同作为活性基因增强子的标志，在边界附近的基因区间内有较为平稳的富集状态；而H3K27me3 和H3K9me3 与基因抑制有关，它们在边界处的峰值变化不明显，甚至有低谷（或损耗）出现。

图5 接触域边界附近CTCF和8种组蛋白峰值富集状态Figure 5 The peak enrichment status of CTCF and proteins near the topological domain boundary

尽管大多数边界富集了CTCF结合位点，但单一变量CTCF 还不足以识别域边界。以基因组峰值分布信息为例，人类H1 细胞系中的22 条常染色体上共有CTCF 峰值63 863 个，只有大约30.26%的CTCF 峰值位点位于接触域边界的附近区域，这种相对较低的数据利用率可能与识别接触域时设置的严格阈值有关。然而，CTCF 不存在的基因位点处也可能会有接触域存在，因此需要借助其他协变量来完成截止阈值的设定。在边界处有显著变化（明显增加或损耗）的组蛋白共有6 种，其中CTCF 和H4K20me1 的组合最高，将CTCF的利用率提高到了35.25%，而CTCF和H4K9me3也在一定程度上将利用率提高到了33.05%，说明缺乏H3K9me3也可以作为识别域边界的重要特征（图6）。

图6 CTCF与各组蛋白修饰组合对数据利用率的影响Figure 6 The effect of CTCF and histone modification combinations on data utilization

在40 kb 的人类H1 细胞系中，分别采用不同的协变量组合方式确定截止阈值，进而调用作用域的边界。表1显示，CTCF作为唯一变量识别的接触域边界数量较多，但只有65.28%的边界是保守性的。一般地，保守性越高说明域边界的预测结果越准确。CTCF 与其他组蛋白的结合会提高截止阈值，使得域边界数量减少，但保守性边界比重增加，表明借助组蛋白修饰的确可以有效提升预测结果的准确性。另外，由于染色体的长度固定，所以相邻域边界之间的距离会随着边界数量的减少而增加。通常两个相邻边界之间的距离越小，说明域结构越精细，能捕捉到嵌套TAD 的概率就越大。

表1 组蛋白修饰对域边界检测结果的影响Table 1 The effect of histone modifications on domain boundary detection

利用类似GSEA 的算法确定截止阈值。图7（a）是经峰值检测识别的所有局部极大峰值的降序排列，图7（b）为降序排列后的峰值位点的富集分数，不同的协变量组合方式对应的富集分数的极大值不同。

图7 GSEA算法确定截止阈值Figure 7 The determination of cut⁃off threshold by GSEA

2.5 性能评价指标

2.5.1 不同域边界检测方法的一致性比较

图8 是HiCID 算法与其他3 种接触域检测算法（HiCDB、TopDom 和DomainCaller［22］）的性能比对分析。表2 和表3 分别统计了人类和小鼠不同细胞系采用不同域边界检测算法的一致性检测结果，以40 kb人类IMR90 为例，由于CTCF 与组蛋白H3K9me3 的组合最高将保守性提高到87.64%，所以该组合被应用于HiCID 算法中以筛选域边界。与hESC 细胞系中的检测结果相似，相较于原HiCDB，NE 降噪后的HiCDB 算法（HiCDB+NE）识别的边界数量有所减少，但边界保守性基本保持稳定。HiCID 算法识别的域边界数量介于HiCDB 算法和HiCDB+NE 之间，但保守性边界占边界总数的百分比均高于此两种算法。尽管TopDom 和DomainCaller识别的域边界数量有限，但保守性边界居多，而且针对不同的物种有较为稳定的检测性能，但是这两种算法计算的相邻边界间的平均距离较大，不能识别更加精细的接触域结构。同时发现，HiCID 算法应用于10 kb 的IMR90 细胞系时，与其他3 种算法相比，不仅识别的边界数量最多、保守性最高，而且相邻边界之间的距离最接近。由此推测，HiCID 算法对较高分辨率的Hi⁃C数据更有效。

图8 人类hESC和IMR90细胞系中不同域边界检测方法的一致性比较Figure 8 The comparison of consistency methods for different domains of human hESC and IMR90

表2 人类细胞系不同域边界检测方法的一致性比较Table 2 The comparison of consistency methods for different domains of human hESC and IMR90

表3 小鼠细胞系不同域边界检测算法一致性比较Table 3 Comparison of the consistency of different cell boundary detection algorithms in mouse cell lines

2.5.2 域边界检测结果的准确性

高活性的组蛋白修饰可用来标记作用域的起点和终止位点，即域边界位置。因此，通过统计比较边界附近各组蛋白修饰的数量可进一步判断域边界的准确性。图9 显示，在整个统计区间内，对于4 种不同的组蛋白修饰，表征HiCID 算法的蓝色曲线与表征HiCDB+NE 算法的红色曲线在个别区间内有交错重叠，但蓝色曲线的整体幅值略高于红色曲线的幅值。

为了从数值上体现两种算法在性能上的差异，分别在hESC 细胞系的22 条常染色体和小鼠Cortex 细胞系19 条常染色体上计算了边界附近单位bin 内平均峰值的个数。从图10 可以看出，HiCID 算法检测的5 种组蛋白的平均峰值的个数均大于HiCDB+NE 算法检测的平均峰值个数，说明在识别域边界时，HiCID 算法在整个基因组都普遍有效。

图9 4种组蛋白修饰在hESC域边界附近平均峰值个数Figure 9 Average number of peaks of four histone modifications near the boundary of hESC domain

2.5.3 接触域边界检测结果

图11是HiCDB与HiCID识别拓扑域边界的结果对比图。人类hESC 细胞系chr18 染色体的［67 100 000，71 100 000］基因区间约有154个bin（分辨率为40 kb），Hi⁃C热图中用浅蓝色的点标注了两种算法识别的域边界位置，而深蓝色的点表示在21 种人类细胞系中都存在的保守域边界。从图11 可以直观地看出，改进前后HiCDB 识别的域边界有部分重合，尤其对保守性边界标注基本一致，而HiCID不仅能识别边界轮廓清晰的域结构，对域内出现的嵌套TAD也能有效识别。

3 讨论与结论

当前表观遗传学领域涌现出众多域检测算法，这些方法为检测拓扑域及其边界提供了更多的选择，但大多数检测方法在应用过程中都有一定的局限性，比如：有的算法需要对多个参数进行调整；有的不适用于有多种分辨率的Hi⁃C数据；或者算法的可重复性差、运行时间成本较高、检测结果的准确性有待提高。

图10 hESC和Cortex细胞系各组蛋白平均出现的峰值个数Figure 10 The average number of peaks of each histone on 19 chromosomes in mice Cortex and hESC

图11 HiCDB与HiCID接触域边界检测结果Figure 11 The detection results of HiCDB and HiCID contact domain boundaries

与代表性的HiCDB 算法相比较，本文提出的Hi⁃CID算法利用网络增强技术对原始Hi⁃C数据进行降噪预处理，构建绝缘密度统计量确定边界特征，具有不易漏检和冗余性好的特点。针对HiCDB算法在筛选阈值时只考虑单一协变量CTCF，本文在识别过程中增加了组蛋白修饰信息，进一步提升了算法稳健性。此外，HiCID 适用范围广泛，对于不同分辨率下的Hi⁃C数据，最优窗口参数的选择可以利用TADquality有效约束。

由于现有实验数据的局限性，本文提出的算法对边界元件信息的融合程度还不够，如管家基因、tRNAs、SINE 反转录转座子等DNA 元件的分布信息未能引入算法，因此在处理不同物种或不同细胞系时算法性能还不够稳定。此外，除算法本身的局限性之外，目前关于染色质相互作用域的大小和划分还没有一个标准的定义，在同一个作用域中会包含大量嵌套的TAD。此外，在细胞过程中会发生改变和重组，给有效定位域边界带来了挑战。