一种Hi-C数据拓扑相关结构域划分方法及装置

一种hi-c数据拓扑相关结构域划分方法及装置
技术领域
1.本发明涉及计算生物学技术领域，特别涉及一种hi-c数据拓扑相关结构域划分方法及装置。


背景技术：

2.hi-c是一种高通量的染色质构象捕获技术，hi-c技术能够通过甲醛固定、酶切、环化连接、序列打断、建库测序并将它们映射到参考基因组来捕获全基因组的染色体相互作用。hi-c技术利用高通量的测序技术，可以生成全基因组、大规模的染色体内和染色体间相互作用(接触频率)数据，可以描述基因组内的空间结构。对相互作用的dna片段进行更深的测序，则可以获得更高分辨率的hi-c数据，显示更多、更精细的结构信息。随着其技术不断发展成熟和测序价格的大幅下降，它在文献中逐渐成为常规使用的一种研究基因组的技术。
3.真核生物基因组在细胞核中折叠成一个具有明确空间结构区域。近些年来三维基因组研究方兴未艾，随着对hi-c等基因组测序数据的研究，已经逐步发掘出基因组的层次结构。基因组的层次结构由大到小包括：活性和非活性区室(a/b compartment)，拓扑相关结构域(tad)，和染色质环。
4.tad最近在文献中备受关注，其特征是在染色体的基因组位点在结构域内具有高水平的相互作用而与结构域外的相互作用水平较低，这一结构承担重要生物学功能。tad的平均长度约为1mbp，不同的tad之间又绝缘子进行分隔，且tad具有跨细胞系稳定性。tad的一个重要功能是形成基因调控的独立区域，同时与邻近区域隔离开来，基因组tad的划分对于研究基因调控、基因组相互作用和基因组功能非常有用。tad破坏可能导致严重的疾病，如癌症。因此，针对实验获得的大量hi-c数据进行划分tad是一项非常重要的研究，对tad的研究有益于深入理解疾病。
5.现有技术中tad的鉴定方法复杂多样，但其中基于聚类的方法，如clustertad等，大多使用40kb分辨率或更低分辨率的hi-c数据用于tad的划分，对于hi-c数据分辨率的稳健性和可靠性较低。在超高分辨率下，由于矩阵数据稀疏性，使得样本在特征空间内的聚集性较差，聚类的准确性下降。所以针对高分辨率下的tad划分方法是值得研究的课题。


技术实现要素：

6.本发明提供了一种hi-c数据拓扑相关结构域划分方法及装置，以增强高分辨率hi-c数据的tad划分的准确性，实现多层次的tad划分。
7.为解决上述技术问题，本发明提供了如下技术方案：
8.一方面，本发明提供了一种hi-c数据拓扑相关结构域划分方法，包括：
9.获取hi-c基因组测序原始观察数据和用于将原始观察矩阵进行归一化的向量，利用所述向量对所述原始观察数据预处理，得到归一化的hi-c接触矩阵；
10.基于hi-c数据接触频率和空间距离之间的函数关系，将所述hi-c接触矩阵转换为
基因组各个位点之间的距离矩阵；
11.通过预设的最短路径算法计算基因组位点间的最短距离，对于每个位点保留距离最小的部分数值，得到基因组的空间距离图谱；
12.通过预设的聚类算法对所述基因组的空间距离图谱中的位点进行聚类，得到hi-c数据的拓扑相关结构域划分。
13.进一步地，所述hi-c基因组测序原始观察数据的分辨率不小于50kb。
14.进一步地，利用所述向量对所述原始观察数据预处理，包括：
15.对所述hi-c基因组测序原始观察数据中的值除以所述向量中的相应范数因子。
16.进一步地，所述hi-c数据接触频率和空间距离之间的函数关系为：
[0017][0018]
其中，f表示hi-c数据接触频率，d表示空间距离，α为预设的转换参数。
[0019]
进一步地，所述预设的最短路径算法为floyd warshall算法。
[0020]
进一步地，所述对于每个位点保留距离最小的部分数值为使用黄金分割搜索算法确定保留值的个数；所述使用黄金分割搜索算法确定保留值的个数包括：
[0021]
设置保留值数量为使用的分辨率下的hi-c数据中染色体位点总个数的10％，即v
keep
＝length(d)/10，其中d是距离矩阵，length(d)表示距离矩阵的行数也就是染色体位点个数；然后定义一个单峰目标函数：
[0022][0023]
其中，v
tad
表示描述tad结构的值，v
longerange
表示远程的相互作用；
[0024]
在搜索区间内计算单峰目标函数值，通过黄金分割搜索算法不断缩小单峰目标函数的最值的已知范围，从而找到函数的最大值，得到保留值的个数。
[0025]
进一步地，所述预设的聚类算法为dbscan算法；
[0026]
通过聚类算法对所述基因组的空间距离图谱中的位点进行聚类，包括：
[0027]
计算临近样本之间的距离并统计其分布，结果随距离的增加显示为单峰的分布，依据样本距离分布，选择峰值作为邻域半径eps参数；
[0028]
将hi-c基因组测序数据的分辨率和拓扑相关结构域结构的大小相除，计算出半径内最少点数minpts；
[0029]
基于eps和minpts，通过dbscan算法对所述空间距离图谱中的位点聚类。
[0030]
另一方面，本发明还提供了一种hi-c数据拓扑相关结构域划分装置，包括：
[0031]
数据处理模块，用于获取hi-c基因组测序原始观察数据和用于将原始观察矩阵进行归一化的向量，利用所述向量对所述原始观察数据预处理，得到归一化的hi-c接触矩阵；
[0032]
距离图谱计算模块，用于基于hi-c数据接触频率和空间距离之间的函数关系，将所述hi-c接触矩阵转换为基因组各个位点之间的距离矩阵；通过预设的最短路径算法计算基因组位点间的最短距离，对于每个位点保留距离最小的部分数值，得到基因组的空间距离图谱；
[0033]
聚类模块，用于通过预设的聚类算法对所述基因组的空间距离图谱中的位点进行
聚类，得到hi-c数据的拓扑相关结构域划分。
[0034]
再一方面，本发明还提供了一种电子设备，其包括处理器和存储器；其中，存储器中存储有至少一条指令，所述指令由处理器加载并执行以实现上述方法。
[0035]
又一方面，本发明还提供了一种计算机可读存储介质，所述存储介质中存储有至少一条指令，所述指令由处理器加载并执行以实现上述方法。
[0036]
本发明提供的技术方案带来的有益效果至少包括：
[0037]
本发明提供的技术方案通过聚类方法对hi-c数据进行拓扑相关结构域划分。tad接触数据通过转换函数得到空间距离，又利用最短路径算法将高分辨率下稀疏的相互作用数据，获得描述空间近邻性的空间距离图谱。使用空间距离图谱作为聚类输入，满足tad物理上的空间分布特性，从而可以充分利用稀疏的高分辨率hi-c数据得到更准确的tad划分。而且本发明可以使用50kb及更高分辨率下的hi-c数据进行划分，在50kb至5kb分辨率的hi-c数据的tad划分过程具有稳定性，能够得到在不同观察尺度下的tad的层次划分。
附图说明
[0038]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0039]
图1是本发明实施例提供的hi-c数据拓扑相关结构域划分方法的流程图；
[0040]
图2是本发明实施例提供的hi-c数据拓扑相关结构域划分方法的详细流程图；
[0041]
图3是本发明实施例提供的hi-c数据拓扑相关结构域划分装置的系统框图。
具体实施方式
[0042]
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
[0043]
第一实施例
[0044]
本实施例提供了一种hi-c数据拓扑相关结构域划分方法，该方法可以由电子设备实现，该方法的执行流程如图1所示，包括以下步骤：
[0045]
s1，获取hi-c基因组测序原始观察数据和用于将原始观察矩阵进行归一化的向量，利用所述向量对原始观察数据预处理，得到归一化的hi-c接触矩阵；
[0046]
其中，hi-c基因组测序原始观察数据的分辨率为50kb或更高。所述向量用于对hi-c基因组测序原始观察数据的系统偏差进行处理，从而对原始观察数据进行矩阵平衡，得到标准化的hi-c接触矩阵。具体地，在本实施例中，用于将原始观察矩阵进行归一化的向量可以使用knight-ruiz(kr)或vanilla-coverage(vc)等常用的hi-c数据归一化方法计算得到。其中，
[0047]
kr归一化方法，根据hi-c接触矩阵具有对称非负的特性，基于此特性进行对角线缩放使用共轭梯度迭代方法将矩阵平衡，使其满足矩阵的每一行和每一列元素之和为1；
[0048]
vc归一化方法对hi-c接触矩阵的每一个元素首先除以相应行的总和，然后除以相应列的总和，得到结果为归一化的矩阵。
[0049]
具体地，如图2所示，本实施例利用所述向量对原始观察数据预处理为：对hi-c基因组测序原始观察数据中的值除以kr归一化向量中的相应范数因子，进行去除噪声操作，得到归一化的hi-c接触矩阵f。
[0050]
s2，基于hi-c数据接触频率和空间距离之间的函数关系，将所述hi-c接触矩阵转换为基因组各个位点之间的距离矩阵；
[0051]
其中，hi-c数据接触频率和空间距离之间满足如下的函数关系：
[0052][0053]
其中，f表示hi-c数据接触频率，d表示空间距离，α为转换参数。该转换参数可以根据分辨率大小进行选择，其默认取值为1/3。
[0054]
基于上述函数关系，将接触矩阵f转换为距离矩阵d。
[0055]
s3，通过预设的最短路径算法计算基因组位点间的最短距离，对于每个位点保留距离最小的部分数值，得到基因组的空间距离图谱；
[0056]
如图2所示，最短路径算法为floyd warshall，用于解决带权图的最短路径距离，其通过循环遍历空间距离矩阵中的顶点，对于其中距离值为无穷远的位置检查是否可以通过基点使一对顶点简单的距离更小，则更新两点间的距离。
[0057]
floyd warshall最短路径算法能够计算基因组位点间的最短距离，对于高分辨率hi-c数据，由于其测序深度的增加，部分空间临近的位点间不存在测序读对，所以高分辨率hi-c数据非常稀疏。通过这一计算步骤，得到矩阵可反映局部空间位点间的近邻性，以便后续样本间距离计算。同时应用保留部分最大值的方法，可避免在应用最短距离算法后的远程相互作用对tad的划分造成影响。
[0058]
具体地，每个位点保留距离最小的部分数值的方法是使用黄金分割搜索算法确定保留值的个数。其中，需要说明的是，对于不同分辨率下需要保留不同数量的点，以满足能够描述tad的结构信息，同时避免远程相互作用对于tad边界划分的影响。基于此，设置保留值数量为使用的分辨率下的hi-c数据中染色体位点总个数的10％，即v
keep
＝length(d)/10，其中，d是距离矩阵，length(d)表示距离矩阵的行数也就是染色体位点个数。然后定义一个单峰目标函数：
[0059][0060]
其中，v
tad
表示可以描述tad结构的值，两位点间的线性距离相聚5mbp以内。v
longerange
表示远程的相互作用，两位点间的线性距离相聚10mbp以上。
[0061]
在搜索区间内计算函数值，若σ未取得最大值，则通过黄金分割搜索算法迭代缩小单峰函数的最值的已知范围(搜索范围)，搜索新的保留值个数，并依据此值重新计算σ，直至找到函数的最大值，得到保留值的个数。
[0062]
s4，通过预设的聚类算法对所述基因组的空间距离图谱中的位点进行聚类，得到hi-c数据的拓扑相关结构域划分。
[0063]
其中，如图2所示，聚类算法为基于密度的聚类方法dbscan；其中dbscan算法的参数选择使用空间距离图谱作为输入其中邻域半径eps参数，使用黄金分割算法，依据给定范
围进行搜索，查找最优参数；半径内最少点数minpts依据所述hi-c基因组测序数据的分辨率和tad结构的大小，相除计算得到。具体地，通过聚类算法对空间距离图谱中的位点进行聚类，包括：
[0064]
计算临近样本之间的距离并统计其分布，结果随距离的增加显示为单峰的分布，依据样本距离分布，选择峰值作为邻域半径eps参数；则样本位点p的eps邻域定义为：n
eps
(p)＝{q∈d|dist(p,q)≤eps}；
[0065]
其中，在本实施例中，基于密度的聚类方法dbscan的样本距离计算函数使用夹角余弦距离dist(p,q)＝cosθ。在二维空间中，向量a(x1,y1)和b(x2,y2)的夹角余弦公式为将其扩展到n维空间，则对于两个n维样本点p(x
11
,x
12
,
…
,x
1n
)和q(x
11
,x
12
,
…
,x
1n
)来说，其计算公式为：
[0066][0067]
依据hi-c基因组测序数据的分辨率和拓扑相关结构域结构的大小，相除计算出半径内最少点数minpts；例如选用数据分辨率为5kb，欲查看染色质长度为500kb或以上大小的拓扑相关结构域时，所选择的minpts应为100。如果给定对象eps邻域内的样本点数大于等于minpts，则称该对象为核心对象。
[0068]
通过合理设置dbscan算法的上述两个重要参数，对所述基因组的空间距离图谱中的位点聚类，得到分类的标签并展示。具体过程如下：
[0069]
步骤1，检查样本点p，如果p未被处理，则检查其邻域，如果样本点p的邻域中包含点数不小于半径内最少点数minpts，则建立新的簇c，邻域半径中的所用点加入候选集n。如果样本点p的邻域中包含点数小于半径内最少点数minpts，则将这个样本标记为噪声点。
[0070]
步骤2，对候选集n中所有尚未被处理的对象q，检查其邻域，如果其中包含点数不小于minpts，则将这些对象加入n，如果q未归入任何一个簇，则将q加入c簇中。
[0071]
步骤3，重复上述步骤2，继续检查n中未被处理的对象，直至候选集n为空，则停止。
[0072]
重复上述步骤1至步骤3，当所有样本都归入某个簇或标记为噪声，则停止。
[0073]
综上，本实施例的技术方案通过聚类方法对hi-c数据进行拓扑相关结构域划分。tad接触数据通过转换函数得到空间距离，又利用最短路径算法将高分辨率下稀疏的相互作用数据，获得描述空间近邻性的空间距离图谱。使用空间距离图谱作为聚类输入，满足tad物理上的空间分布特性，从而可以充分利用稀疏的高分辨率hi-c数据得到更准确的tad划分。而且本发明可以使用50kb及更高分辨率下的hi-c数据进行划分，在50kb至5kb分辨率的hi-c数据的tad划分过程具有稳定性，能够得到在不同观察尺度下的tad的层次划分。
[0074]
第二实施例
[0075]
本实施例提供了一种hi-c数据拓扑相关结构域划分装置，该hi-c数据拓扑相关结构域划分装置的系统结构如图3所示，包括以下模块：
[0076]
数据处理模块，用于获取hi-c基因组测序原始观察数据和用于将原始观察矩阵进行归一化的向量，利用所述向量对所述原始观察数据预处理，得到归一化的hi-c接触矩阵；
[0077]
距离图谱计算模块，用于基于hi-c数据接触频率和空间距离之间的函数关系，将所述hi-c接触矩阵转换为基因组各个位点之间的距离矩阵；通过预设的最短路径算法计算
基因组位点间的最短距离，对于每个位点保留距离最小的部分数值，得到基因组的空间距离图谱；
[0078]
聚类模块，用于通过预设的聚类算法对所述基因组的空间距离图谱中的位点进行聚类，得到hi-c数据的拓扑相关结构域划分。
[0079]
本实施例的hi-c数据拓扑相关结构域划分装置与上述第一实施例的hi-c数据拓扑相关结构域划分方法相对应；其中，本实施例的hi-c数据拓扑相关结构域划分装置中的各功能模块所实现的功能与上述第一实施例的hi-c数据拓扑相关结构域划分方法中的各流程步骤一一对应；故，在此不再赘述。
[0080]
第三实施例
[0081]
本实施例提供一种电子设备，其包括处理器和存储器；其中，存储器中存储有至少一条指令，所述指令由处理器加载并执行，以实现第一实施例的方法。
[0082]
该电子设备可因配置或性能不同而产生比较大的差异，可以包括一个或一个以上处理器(central processing units，cpu)和一个或一个以上的存储器，其中，存储器中存储有至少一条指令，所述指令由处理器加载并执行上述方法。
[0083]
第四实施例
[0084]
本实施例提供一种计算机可读存储介质，该存储介质中存储有至少一条指令，所述指令由处理器加载并执行，以实现上述第一实施例的方法。其中，该计算机可读存储介质可以是rom、随机存取存储器、cd-rom、磁带、软盘和光数据存储设备等。其内存储的指令可由终端中的处理器加载并执行上述方法。
[0085]
此外，需要说明的是，本发明可提供为方法、装置或计算机程序产品。因此，本发明实施例可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本发明实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质上实施的计算机程序产品的形式。
[0086]
本发明实施例是参照根据本发明实施例的方法、终端设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、嵌入式处理机或其他可编程数据处理终端设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理终端设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0087]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理终端设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。这些计算机程序指令也可装载到计算机或其他可编程数据处理终端设备上，使得在计算机或其他可编程终端设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程终端设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0088]
还需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要
素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。
[0089]
最后需要说明的是，以上所述是本发明优选实施方式，应当指出，尽管已描述了本发明优选实施例，但对于本技术领域的技术人员来说，一旦得知了本发明的基本创造性概念，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明实施例范围的所有变更和修改。