基于谱聚类的染色质拓扑关联结构域预测方法及电子装置

1.本发明涉及生物信息技术领域，具体涉及一种基于谱聚类的染色质拓扑关联结构域预测方法及电子装置。


背景技术：

2.dna的复制、基因调控和基因表达等生物学功能都依赖于染色质三维结构实现。研究表明，特定基因组区域和特定条件下染色体结构的变化与多种人类疾病高度相关，包括癌症、发育过程中肢体畸形和严重的脑部发育异常。因此，研究染色质三维结构对于解释基因的表达和调控等细胞过程，进而指导遗传病相关研究和基因治疗等医学问题方面具有重要价值。
3.拓扑关联结构域(topologically associated domain，tad)作为染色质层次结构中的一种，于2012年被首次发现。该区域富集了绝缘子结合蛋白ctcf、管家基因、sine逆转座子等，在基因调控中发挥关键作用。进一步的研究发现，在苍蝇、蠕虫、真菌和细菌中也检测到类似的结构域。因此，这些区域在细胞分裂中是稳定存在的，在不同细胞系中存在一定的进化保守性。tad目前已被认为是染色体折叠的基本单位，并被认为是染色体组织中的一个重要二级结构。
4.然而，传统生物实验识别tad区域存在耗时长、花费成本高、实验实施困难等问题。为更好地理解tad在生物生长发育和遗传过程中发挥的功能，研究者需要可靠的预测模型指导生物实验定位染色质上的tad区域。
5.根据预测模型所用数据的不同，可将该领域已有的计算方法分为两类：基于3c-based数据的预测方法与基于组蛋白修饰信号数据的预测方法。
6.基于3c-based数据的预测方法是最早定义和研究染色质拓扑关联结构域的方法，主要利用全基因组范围上的染色质读取段间的交互频次矩阵。依据各读取段的交互频次的差异，对各读取段进行分类，进而获取tad区域。但本技术发明人在实施本发明的过程中，发现现有基于3c-based数据的染色体拓扑关联结构域计算方法存在一些不足之处：存在假阳性高、参数选择困难等问题。这类方法的实施可以划分为特征提取、模型构建和结果筛选三个阶段。当前特征提取阶段未能很好提取对应的特征向量来量化染色质读取段交互频次之间的差异。模型阶段由于区域划分的错误，导致结果假阳性过高。结果筛选阶段则会因为未能有效剔除假阳性区域，故而难以准确区分tad区域和“gap”区域。
7.基于组蛋白修饰信号数据的方法是利用拓扑关联结构域边界处存在ctcf位点、组蛋白和一些基因调控元件富集或缺失的现象，通过构建计算模型预测拓扑关联结构域的边界。但这类方法存在特征选择困难的问题，难以准确提取有效的特征。


技术实现要素：

8.针对上述问题，本发明的目的在于提供一种基于谱聚类的染色质拓扑关联结构域预测方法及电子装置，从位点(染色质读取段)上下游的角度聚焦位点的交互特征，并使用
余弦相似度量化位点交互模式的相似性，进而提高模型预测准确性。技术方案如下：
9.一种基于谱聚类的染色质拓扑关联结构域预测方法，包括如下步骤：
10.s1：获取人类常见细胞系的hi-c数据，并进行数据预处理；
11.s2：对于hi-c数据中的每个位点，分别提取上下游交互频次数据作为该位点的特征向量；
12.s3：根据提取的特征向量，使用余弦相似度计算位点与位点之间的相似性，构建对应的相似性矩阵；
13.s4：基于相似性矩阵完成相似性无向图的构建，再使用聚类算法对位点进行分类；
14.s5：从聚类结果中提取tad区域并进行筛选，预测tad区域。
15.进一步的，所述s1具体包括：
16.s11：获取hi-c数据并对该数据进行规范化以消除位点之间距离带来的数据噪声；
17.s12:对输入的hi-c频次矩阵进行ln(x+1)处理，以减小hi-c交互频次数据的动态范围，并进行数据平滑，使交互数据更符合高斯分布；并对每个频次加1，以避免出现负无穷大的值。
18.更进一步的，所述步骤s2具体为：
19.s21：对于每个位点分别提取上下游2mb范围的染色质交互频次数据作为该位点的特征向量；对于没有足够的上/下游的位点，用上/下游的所有交互频次数据的平均值进行填充；
20.s22：将得到的特征向量拼接为特征矩阵。
21.更进一步的，所述步骤s3具体为：
22.s31：根据任意两位点的特征向量，计算位点之间的余弦相似性，越大的余弦值表明两个向量方向越接近；余弦相似性计算如下：
[0023][0024]
其中，cosine
ij
表示位点i和位点j的特征向量余弦值；fi和fj分别表示位点i和位点j的特征向量；
[0025]
s32：通过对余弦值进行数值变换，将位点之间的相似性范围变为[0,1]，得到相似性矩阵；相似性矩阵计算公式为：
[0026][0027]
其中，s
ij
为相似性矩阵s中的元素。
[0028]
更进一步的，所述步骤s4具体为：
[0029]
s41：定义为无向图，其中顶点集v＝{v1,v2,
…
,vn}表示位点，边集e表示位点与位点之间的相似性；图中的邻接关系表示对称的相似性矩阵中的邻接关系；
[0030]
s42：由于顶点的度指与该顶点相关联的边的条数，则定义顶点vi的度d(vi)为：
[0031]
[0032]
其中，表示相似度矩阵，顶点vi和顶点vj之间的边的权重，也就是两个顶点之间相似性；
[0033]
度矩阵为对角矩阵，定义为：
[0034][0035]
s43：计算拉普拉斯矩阵，计算公式为：
[0036][0037]
对拉普拉斯矩阵进行规范化：
[0038][0039]
计算规范化的拉普拉斯矩阵中前k个特征值对应的特征向量u1,
…
,uk，将其组成矩阵u∈rn×k，u＝[u1,
…
,uk]；
[0040]
s44：在矩阵u上使用k-means算法，将顶点聚集到c1,
…
,ck，k个簇中。
[0041]
更进一步的，所述步骤s5具体为：
[0042]
s51：定义tad区域：tad的最小尺寸为180kb，将矩阵对角线上同一簇的连续位点连接成段，将属于同一类别且大于180kb的片段预测为tad，将小于180kb的片段定义为tad边界；
[0043]
s52：剔除假阳性区域：若一个tad的域内平均交互频次低于整个染色体的平均交互频次，则该区域为被聚集在一起的“gap”区域，而非拓扑关联结构域区域，将这其从结果中剔除。
[0044]
一种基于谱聚类的染色质拓扑关联结构域预测的电子装置，包括处理器，存储器以及储存在存储器上的计算机程序，所述计算机程序在处理器上执行实现上述方法。
[0045]
本发明的有益效果是：
[0046]
1)本发明将位点上下游2mb范围内的交互频次数据作为该位点的染色质相互作用模式特征，为后续模型提供更多有效的交互特征；
[0047]
2)本发明引入余弦相似度量化位点之间相互作用模式的相似性，相似性较高的两位点之间的边具有更高的权重，在聚类模型中具有更高的区分度，进而提高模型预测准确性；
[0048]
3)在相同数据集上与现有模型对比，取得了最优的预测结果，并在相关生物学数据进行验证，以证明模型预测结果具有生物学意义；
[0049]
4)本发明通过已有hi-c数据，来预测未被验证的区域是否为tad区域以指导生物实验，有效地减少了实验时间与财力损耗。
附图说明
[0050]
图1是本发明基于谱聚类的染色质拓扑关联结构域预测方法和系统的流程图。
[0051]
图2是本发明提出的方法与该领域的三种典型方法在(a)tad质量、(b)皮尔森相关系数和(c)距离相关系数三种指标上的性能对比图。
[0052]
图3是本发明提出的方法与该领域的三种典型方法在不同噪声水平的数据上在(a)皮尔森相关系数和(b)距离相关系数两种指标上的性能对比图。
[0053]
图4是本发明提出的方法预测tad边界组蛋白修饰信号富集程度分析结果图；其中，(a)-(f)属于imr90人类细胞系，(g)-(l)属于hesc人类细胞系。
[0054]
图5是本发明提出的方法预测tad边界的motif富集程度分析结果图；(a)总目标序列＝7601，总背景序列＝7681；(b)总目标序列＝6370，总背景序列＝6414。
具体实施方式
[0055]
下面结合附图和具体实施例对本发明做进一步详细说明。
[0056]
本发明提出了一种基于谱聚类的染色质拓扑关联结构域预测方法和系统。该方法通过提取位点上下游2mb范围内的交互频次数据表示位点上下文的交互频次特征关系，并引入余弦相似度量化位点之间交互模式的相似性，赋予相似度高的两位点之间的边更大的权重，使用谱聚类方法构建预测模型，进而提高模型的预测准确性。
[0057]
本实施例提供了一种基于谱聚类的染色质拓扑关联结构域预测方法，参考图1与图2，其过程是基于python3.8.6-scikit-learn0.24.0实现。该方法包括：
[0058]
s1：数据获取与预处理。获取人类细胞系常见的hi-c数据，并进行预处理；
[0059]
s2：位点特征向量提取。分别提取位点上下游交互频次数据作为位点特征向量；
[0060]
s3：相似度矩阵的构建。根据提取后的特征向量，使用余弦相似度计算位点与位点之间的相似性，构建相似度矩阵；
[0061]
s4：谱聚类。依据相似度矩阵，使用聚类算法对位点进行分类；
[0062]
s5：tad区域的识别。从聚类结果中识别tad区域并进行筛选。
[0063]
具体而言，上述五个步骤的详细过程为：
[0064]
一、数据获取与预处理
[0065]
从公开数据库中获取hi-c数据集，并进行规范化和预处理。
[0066]
1、获取数据。获取的hi-c数据来源于美国国家生物信息中心(ncbi)公开数据库geo中的数据集，访问号为gse35156，数据集的分辨率为40kb。使用yaffe等人提供的规范化方法处理数据，以消除生物实验中由于位点之间间隔的距离因素产生的偏差。
[0067]
2、数据预处理。对输入的hi-c频次矩阵进行ln(x+1)处理，以减小hi-c交互频次数据的动态范围，并进行数据平滑，使hi-c数据更符合高斯分布。由于两个位点之间的交互频次可能为零，因此需要对每个频次加1，以避免出现负无穷大的值。
[0068]
二、特征向量提取
[0069]
1、对于预处理好后的数据，提取每个位点上下游各2mb范围内的相互作用频次数据。数据集的分辨率为40kb，等效于提取位点与上下游各50个位点之间的交互频次数据，最终得到的特征向量长度为101。对于起始的50个位点/最后的50个位点，其没有足够的上/下游位点，需要使用数据填充策略。本发明用上/下游的所有交互频次数据的平均值进行填充，以固定所有位点的特征向量长度均为101。
[0070]
三、相似度矩阵构建
[0071]
使用提取的位点的特征向量计算位点与位点之间的相似度。本发明使用余弦相似度量化位点的相似度，计算如下：
[0072][0073]
余弦相似度更关注数据之间的模式而不是数量的大小，非常适用于量化位点交互模式的相似性。余弦值的范围为[-1,1]，越大的余弦值表明两个向量方向越接近。通过对余弦值进行变换，将位点之间的相似性范围变为[0,1]，得到相似性矩阵s，越大的相似性值，意味着两个位点之间的交互模式越相似。范围变换计算公式：
[0074][0075]
四、聚类模型
[0076]
构建完成相似性矩阵后，使用谱聚类算法对位点进行分类。
[0077]
定义为无向图，其中顶点集v＝{v1,v2,
…
,vn}表示位点，边集e表示位点与位点之间的相似性。图中的邻接关系表示对称的相似性矩阵s中的邻接关系。谱聚类的计算过程如下：
[0078]
一个顶点的度是指与该顶点相关联的边的条数。给的顶点vi的度d(vi)定义为：
[0079][0080]
其中，表示相似度矩阵，顶点vi和顶点vj之间的边的权重，也即两个顶点之间相似性；
[0081]
度矩阵为对角矩阵，定义为：
[0082][0083]
拉普拉斯矩阵定义为：
[0084][0085]
规范化的拉普拉斯矩阵：
[0086][0087]
计算规范化的拉普拉斯矩阵的前k个特征值对应的特征向量u1,
…
,uk，将其组成矩阵u∈rn×k，u＝[u1,
…
,uk]。在矩阵u上使用k-means算法，将顶点划分到c1,
…
,ck，共k个簇中。
[0088]
五、tad区域识别和筛选
[0089]
1.tad区域识别。tad被视为是有相互交互模式的位点的集群。因此，本发明从上一步获取的聚类簇中识别tad区域。首先将矩阵对角线上属于同一类别的连续位点连接成段，根据以往的文献报道，tad的最小尺寸为180kb，因此将属于同一类别并且大于180kb的片段预测为tad，小于180kb的片段被定义为tad边界。
[0090]
2.tad区域筛选。基于tad的性质，tad内部的染色质平均交互作用频次一般高于相邻两个tad之间的平均交互频次。如果一个tad的域内平均交互频次低于整个染色体的平均
交互频次，那么表明该区域应该是被聚集在一起的“gap”区域，而非tad区域。因此，本发明从识别结果中剔除“gap”区域，以消除假阳性区域。
[0091]
本发明提供的预测方法，在具体实施时，可采用软件方式实现流程的自动运行。运行流程的装置也应当在本发明的保护范围内。
[0092]
以下通过对比实验来验证本发明的有益效果。
[0093]
本实验采用的数据从公开数据库geo中提取而得，包括两类人类细胞系(人胚胎干细胞hesc和人胚肺成纤维细胞imr90)和两类小鼠细胞系(小鼠大脑皮层细胞mcortex和小鼠胚胎干细胞mesc)的hi-c数据集。此外，还包括一个包含不同噪声水平的模拟hi-c数据集，用于验证本发明的方法在不同噪声情况下的稳健性。
[0094]
我们将本发明提出的方法sbtd与该领域的三种典型方法进行了比较，它们分别是topdom、di和spectraltad。评估指标使用tad质量、皮尔森相关性系数(pearson correlation coefficient)和距离相关性系数(distance correlation)。其中，tad质量又称为域内域间交互频次差异(inter-intra diff)，计算公式：
[0095]
diffi＝intra(i)-inter(i,j),|i-j|＝1
[0096]
距离相关性系数计算公式：
[0097][0098][0099][0100][0101]
x
ij
＝||x
i-xj||；y
ij
＝||y
i-yj||
[0102]
上述三种指标，均是数值越高代表方法性能越优秀。从图2可见，本发明方法的结果不管在最大值、中位数值还是最小值均优于上述三种方法，表明本发明的方法预测的tad区域具有更好的内聚性，更接近物理上真实的tad区域。此外，我们还在具有不同噪声水平的hi-c数据集上将本发明的方法和上述三种方法进行了稳健性比较。从图3可见，随着噪声水平的增加，四种方法的性能均有所下降，但是本发明的方法均优于上述三种方法。
[0103]
除了对比方法性能外，我们还验证了本发明方法预测的tad区域是否具有生物学意义。图4为本发明方法预测的tad区域边界的组蛋白修饰信号富集程度结果。tad边界处有特定组蛋白修饰信号富集或缺失的特征现象。图4显示的结果与现有文献总结的结果类似，这一定程度证明了本发明方法的有效性。
[0104]
更进一步的，我们使用homer软件对本发明方法识别的tad边界(target sequences)和非边界(background sequences)进行了motif富集分析。motif是指dna序列中局部的保守区域，是一组序列中共有的一小段序列模式，具有特定的生物学调控功能。从图5可见，ctcf、boris、sp1、sox17、hinfp、e2f4、pgr motifs在tad边界中显着富集。ctcf介
导染色质环的形成并参与tad边界的拓扑约束。boris是ctcf的旁系同源物之一，归类为绝缘子蛋白，在染色质组织和基因表达的调控中起重要作用。其余motif或多或少和基因的抑制表达相关。这与tad边界具有抑制域间基因调控的功能符合。
[0105]
由此可得出结论，与已有染色质拓扑关联结构域预测方法相比，本发明方法拥有更好的性能。
[0106]
综上所述，本发明设计了一种基于谱聚类的染色质拓扑关联结构域预测方法和系统，能够有效提高预测tad区域的性能。本发明的研究成果可应用于生物医学领域，研究人员可以通过本方法的帮助，更好地识别tad区域，进而对tad在基因调控、疾病产生机理和精准医疗等方面进行更深入地研究。此外，由于余弦相似度能更好的量化交互模式的相似性，本发明研究成果不仅能应用于染色质拓扑关联结构域预测问题中，还能应用于其他发生相互作用的基因元件之间的预测。