CN108197431B - 染色质相互作用差异的分析方法和系统 - Google Patents

Abstract

本发明涉及一种不同状态的样品中染色质相互作用差异的分析方法，更具体涉及一种结合近邻相互作用的信息进行差异性染色质相互作用分析的方法。

Description

染色质相互作用差异的分析方法和系统

技术领域

本发明涉及一种不同状态的样品中染色质相互作用差异的分析方法，更具体涉及一种结合近邻相互作用的信息进行差异性染色质相互作用分析的方法。

背景技术

伴随着人类基因组计划的完成，染色质一维结构的信息即DNA序列信息逐步被破解与完善。大量开展的全基因组或外显子组测序等工作，为基因组的序列信息、SNV、功能元件等基因组的研究奠定基础。进一步的，ATAC-seq、ChIP-exo、ChIRP-seq等技术的出现则能够有效的揭示基因组开放程度、转录因子结合情况等信息，推动研究者构建起包含基因-基因、基因-蛋白质相互作用在内的二维调控网络。基于近端连接技术的3C、4C、5C、Hi-C、ChIA-PET和DamID等技术的出现，揭示了蛋白质介导的染色质长程相互作用信息，打开了3D基因组结构研究的大门。而single cell Hi-C、Dilution Hi-C、in situ Hi-C、DNase Hi-C、Capture-C等技术的演变，又使得三维基因组研究的信噪比和精度不断提。

在很多生物学过程中，例如细胞分化、病变，染色质三维结构或者构象均发生了巨大的改变，这就需要对不同时间下和/或状态下所出现的差异染色质相互作用进行分析。但是目前染色质相互作用信号的分析方法，从准确度到精度都有所不足，例如最简单也是最常用的策略之一是使用染色质相互作用频率的变化倍数作为检测的标准，这常常会导致分析结果出现较高数量的假阳性或假阴性。因此本领域迫切需要一种新的方法，能够更加有效的识别染色质相互作用的差异。

发明内容

发明人经过长期的研究，发明了一个染色质相互作用差异分析的新方法,该方法基于线性距离相隔较远的两个染色质区段位置在空间上发生相互作用(或者相互靠近)时，这个两个DNA位置各自周边的染色质区段位置也会出现相互靠近的因素，将周边染色质区段的相互接触强度也纳入了分析范围，并通过合理的模型构建从而识别染色质相互作用的差异。结果发现本发明的方法能够更好的符合染色质相互作用的特点。通过模拟和实际实验数据的分析表明,本发明方法相对目前广泛应用的count-based算法(假设相互作用的DNA位置与其附近的其他DNA位置是相对孤立的),具有更好的检出效果。

本发明至少包括下列实施方案：

1、一种用于分析不同状态的样品中染色质相互作用差异的方法，其包括下列步骤：

(1)分别对处于第一状态和第二状态的样品进行染色质相互作用的分析，

(2)根据作用染色质区段对在第一状态和第二状态的样品间的接触强度差异以及邻居染色质区段对之间的接触强度差异，分析得到在不同状态的样品中具有差异的染色质相互作用，

所述作用染色质区段对是指经过染色质相互作用分析被鉴定实际发生染色质相互作用的染色质区段对，其包括至少两个作用染色质区段；

所述邻居染色质区段对是由邻居染色质区段两两组合形成的，所述邻居染色质区段是指发生相互作用的两个染色质区段各自所临近的染色质区段，由于作用染色质区段对的关系，作用染色质区段附近的染色质区段之间的空间距离也被拉近而产生染色质接触的信号。

优选的，步骤(1)的染色质相互作用分析重复进行n_c次实验，其中c＝1或2，n₁表示针对第一状态的样品重复的次数，n₂表示针对第二状态的样品重复的次数，n₁和n₂相等或不相等，n_c是1-100之间的自然数，例如1、2、3、4、5、6、7、8、9、10、11、12、13、14、15、16、17、18、19或20，优选2、3、4或5。

其中，优选的，步骤(1)的分析能够获得全基因组范围内或者目标区域范围内所存在的一个或更多个染色质相互作用，

更优选的，步骤(1)的分析能够获得全基因组范围内或者目标区域范围内所存在的染色质相互作用位置和强度，

2、根据实施方案1所述的方法，其中所述步骤(2)中获得染色质相互作用的差异包括下列步骤：

a)计算处于第一状态的样品中，每一个或感兴趣的作用染色质区段对在n₁次重复实验中的平均染色质作用强度；和/或处于第二状态的样品中，每一个或感兴趣的作用染色质区段对在n₂次重复实验中的平均染色质作用强度；

b)以步骤a)获得的处于第一状态的样品或处于第二状态的样品的平均染色质作用强度为基准强度，并为该基准强度赋予与所述作用染色质区段相同的空间位置信息，分别计算每一次重复实验中，处于第一状态和第二状态的样品中每一个或者感兴趣的作用染色质区段对及其邻居染色质区段对相对于基准强度的排序，所述排序综合了相对空间位置和相对染色质接触强度的信息，优选的，所述排序是基于各染色质区段对相对于基准强度的欧式距离而确定的，

c)分别获得针对每个样品在不同次的实验中具有相同排序的染色质区段对的集合，建立上述集合相对于基准强度的分布，分析每一排序下的所述分布在处于第一状态的样品和第二状态的样品之间是否具有显著性差异，

优选的，步骤c)中的分布是利用数学模型计算获得，更优选的，所述数学模型是空间泊松点过程(空间泊松分布)，更优选的，基于该空间泊松点过程计算得到该分布下观察到具有该排序的染色质接触强度的概率，并分析上述概率值在第一状态的样品和第二状态的样品之间是否具有显著性差异。

3、根据实施方案2所述的方法，其中分别以处于第一状态和第二状态的样品的平均染色质作用强度为基准强度，各进行一次空间泊松点过程的计算。

4、根据实施方案2或3所述的方法，其中分析所述分布之间是否具有显著性差异时采用Fisher检验，获得p值，随后将由各个排序下的分布获得的p值集合进行rOP统计(r-thordered p-value statistics)，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

5、根据实施方案1至4中任一项所述的方法，所述邻居染色质区段对由位于作用染色质区段两侧各1-10个位置的邻居染色质区段之间两两组合形成，优选1个位置、2个位置或3个位置，更优选1个或2个位置。

6、根据实施方案1至5中任一项所述的方法，其中步骤(2)按照下述步骤进行：

i)为一个作用染色质区段对赋予坐标信息(i,j)，其中i表示发生相互作用的两个染色质区段中，位于基因组序列上游的染色质区段位置，j表示位于基因组序列下游的染色质区段位置；以(i,j)为平面中心建立一个二维平面窗口，宽度为W，所述窗口共包含W²个染色质区段对，用w表示位于与i或j的坐标距离，则上述窗口内各染色质区段对的横坐标可表示为(i+w)、i、和(i-w)，纵坐标为(j+w)、j和(j-w)，除了中心位置(i,j)外，其他位置作为邻居位置，表示为(i[w],j[w])；其中w<＝(W-1)/2，W为3-19之间的奇数，例如3、5、7、9、11、13、15、17、19，优选3、5或7。

ii)将第m次重复试验中(i,j)位置的染色质相互作用频率记做Z(i，j)_m，其中m是自然数，并且m<＝n_c；将位置坐标(i,j)以及该位置的频率值转换为一个三维坐标系的点(i，j，Z(i，j)_m)；同时将窗口内其他位置的染色质接触强度记做Z(i[w],j[w])_m，进一步结合位置坐标将染色质区段对转换为多个三维坐标系的点(i[w]，j[w]，Z(i[w],j[w])_m)；

iii)计算在处于第一种状态的样品中(i，j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽¹⁾)，并在前述三维坐标系中指定一个新的点值μ₁,其坐标为[i，j，mean(Z(i，j)_n ⁽¹⁾)]；和/或

计算处于第二种状态下的样品中(i,j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽²⁾)，并在前述三维坐标系中指定一个新的点值μ₂，其坐标为[i，j，mean(Z(i，j)_n ⁽²⁾)]

iv)分别计算处于第一种状态和第二次状态下的样品的第m次重复实验中，包含中心位置(i，j，Z(i，j)_m)以及邻居位置(i[w]，j[w]，Z(i[w],j[w])_m)在内的点到μ₁或μ₂的欧式空间距离，从而获得第一状态和第二状态下分别与点μ₁或μ₂的欧式空间距离最近的点，称作第一邻居点，空间距离第二近的点称作第二邻居点，以此类推，空间距离第W²近的点称作第W²邻居点，将第一状态下第m次重复实验中的第k邻居点指定为P_mk ⁽¹⁾、将第二状态下第m次重复实验中的第k邻居点被记做P_mk ⁽²⁾，其中k<＝W²；

v)基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k，从而获得W²个判定标准p_k；和/或

基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中μ₂到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k’，从而获得W²个判定标准p_m’；

vi)利用rOP方法统计上述(W)²个p_m值和/或(W)²个p_m’值，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

7、根据实施方案6所述的方法，其中所述的空间泊松分布模型的分析如下：

将点P_mk ⁽¹⁾与P_mk ⁽²⁾的接触强度分别使用

和

表示，其中

如果比值

明显偏离整体分布，可认为在μ₁处第k邻居点的强度发生变化。在零假设

的情况下，比率R_k(μ₁)分别遵循具有2n₁k和2n₂k自由度的Fisher分布。

8、根据实施方案6或7所述的方法，其中所述rOP检验基于以下假设：

其中，θ_k是检验k的有效大小.S_r是p值的第r阶统计量，并且

S_r～Beta(r，2W²-r+1)，其中r值为小于p*W²的最大整数值，其中p∈(0,1]，具体含义表示两个条件之间一个染色质相互作用被视为差异性染色质相互作用时，所含有的产生显著变化的第k邻居的百分比。

9、一种鉴定调控染色质相互作用的试剂的方法，其包括将使样本与一种或多种试剂接触，利用实施方案1至8中任一项所述的方法分析染色质相互作用差异，以及

鉴定相比于不添加调控试剂的对照组能够改变相互作用的试剂。

10、一种分析细胞分化、发育或病变过程中遗传物质高级结构改变的方法，其包括实施方案1-8任一项所述的步骤。

11、一种鉴定染色质结构变异的方法，其包括实施方案1至8中任一项所述的步骤。

12、一种鉴定能够调控遗传物质高级结构或引起染色质结构变异的调控试剂的方法，其包括将使样本与一种或多种试剂接触，利用实施方案1至8中任一项所述的方法分析染色质相互作用差异，以及

鉴定相比于不添加调控试剂的对照组能够改变相互作用的试剂。

13、一种用于鉴定一种或多种指示特定疾病状态的核苷酸相互作用的方法，其包括进行实施方案1至8任一项的步骤，其中在步骤(1)中，提供病人和健康样品，显示有差异的染色质相互作用指示所述相互作用可用于指示特定的疾病状态；所述疾病优选是遗传疾病、免疫性疾病或癌症。

14、一种与染色质结构改变相关的疾病的诊断方法，其包括进行实施方案1至8中任一项的步骤，其中步骤(1)包括提供来自受试者的样品，并根据染色质相互作用的结果判断是否可能患有疾病；所述疾病优选是遗传疾病或癌症。

15、一种染色质相互作用差异的鉴定系统，所述系统包括染色质相互作用实验模块和数据分析模块，其中所述染色质相互作用实验模块用于分别对处于第一状态和第二状态的样品中所存在染色质相互作用进行分析；

所述数据分析模块用于根据作用染色质区段对在第一状态和第二状态的样品间的接触强度差异以及邻居染色质区段对之间的接触强度差异，分析得到在不同状态的样品中具有差异的染色质相互作用；

所述作用染色质区段对是指经过染色质相互作用分析被鉴定实际发生染色质相互作用的染色质区段对，其包括至少两个作用染色质区段；

所述邻居染色质区段对是由邻居染色质区段两两组合形成的，所述邻居染色质区段是指发生相互作用的两个染色质区段各自所临近的染色质区段，由于作用染色质区段对的关系，作用染色质区段附近的染色质区段之间的空间距离也被拉近而产生染色质接触的信号。

优选的，其中染色质相互作用的分析重复进行n_c次实验，其中c＝1或2，n₁表示针对第一状态的样品重复的次数，n₂表示针对第二状态的样品重复的次数，n₁和n₂相等或不相等，n_c是1-100之间的自然数，例如1、2、3、4、5、6、7、8、9、10、11、12、13、14、15、16、17、18、19或20，优选2、3、4或5。

16、根据实施方案15所述的系统，其中数据分析模块获得染色质相互作用的差异包括下列步骤：

a)计算处于第一状态的样品中，每一个或感兴趣的作用染色质区段对在n₁次重复实验中的平均染色质作用强度；和/或处于第二状态的样品中，每一个或感兴趣的作用染色质区段对在n₂次重复实验中的平均染色质作用强度；

b)以步骤a)获得的处于第一状态的样品或处于第二状态的样品的平均染色质作用强度为基准强度，并为该基准强度赋予与所述作用染色质区段相同的空间位置信息，分别计算每一次重复实验中，处于第一状态和第二状态的样品中每一个或者感兴趣的作用染色质区段对及其邻居染色质区段对相对于参比强度的排序，所述排序综合了相对空间位置和相对染色质接触强度的信息，优选的，所述排序是基于各染色质区段对相对于基准强度的欧式距离而确定的，

c)分别获得针对每个样品在不同次的实验中具有相同排序的染色质区段对的集合，建立上述集合相对于基准强度的分布，分析每一排序下的所述分布在处于第一状态的样品和第二状态的样品之间是否具有显著性差异，

优选的，步骤c)中的分布是利用数学模型计算获得，更优选的，所述数学模型是空间泊松点过程(空间泊松分布)，更优选的，基于该空间泊松点过程计算得到该分布下观察到具有该排序的染色质接触强度的概率，并分析上述概率值在第一状态的样品和第二状态的样品之间是否具有显著性差异。

17、根据实施方案16所述的系统，其中分别以处于第一状态和第二状态的样品的平均染色质作用强度为基准强度，各进行一次空间泊松点过程的计算。

18、根据实施方案16或17所述的系统，其中分析所述分布之间是否具有显著性差异时采用Fisher检验，获得p值，随后将由各个排序下的分布获得的p值集合进行rOP统计(r-th ordered p-value statistics)，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

19、根据实施方案15-19中任一项所述的系统，所述邻居染色质区段对由位于作用染色质区段两侧各1-10个位置的邻居染色质区段之间两两组合形成，优选1个位置、2个位置或3个位置，更优选1个或2个位置。

20、根据实施方案15至19中任一项所述的系统，其中步骤(2)按照下述步骤进行：

i)为一个作用染色质区段对赋予坐标信息(i,j)，其中i表示发生相互作用的两个染色质区段中，位于基因组序列上游的染色质区段位置，j表示位于基因组序列下游的染色质区段位置；以(i,j)为平面中心建立一个二维平面窗口，宽度为W，所述窗口共包含W²个染色质区段对，用w表示位于与i或j的坐标距离，则上述窗口内各染色质区段对的横坐标可表示为(i+w)、i、和(i-w)，纵坐标为(j+w)、j和(j-w)，除了中心位置(i,j)外，其他位置作为邻居位置，表示为(i[w],j[w])；其中w<＝(W-1)/2，W为3-19之间的奇数，例如3、5、7、9、11、13、15、17、19，优选3、5或7。

ii)将第m次重复试验中(i,j)位置的染色质相互作用频率记做Z(i，j)_m，其中m是自然数，并且m<＝n_c；将位置坐标(i,j)以及该位置的频率值转换为一个三维坐标系的点(i，j，Z(i，j)_m)；同时将窗口内其他位置的染色质接触强度记做Z(i[w],j[w])_m，进一步结合位置坐标将染色质区段对转换为多个三维坐标系的点(i[w]，j[w]，Z(i[w],j[w])_m)；

iii)计算在处于第一种状态的样品中(i，j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽¹⁾)，并在前述三维坐标系中指定一个新的点值μ₁,其坐标为[i，j，mean(Z(i，j)_n ⁽¹⁾)]；和/或

计算处于第二种状态下的样品中(i,j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽²⁾)，并在前述三维坐标系中指定一个新的点值μ₂，其坐标为[i，j，mean(Z(i，j)_n ⁽²⁾)]

iv)分别计算处于第一种状态和第二次状态下的样品的第m次重复实验中，包含中心位置(i，j，Z(i，j)_m)以及邻居位置(i[w]，j[w]，Z(i[w],j[w])_m)在内的点到μ₁或μ₂的欧式空间距离，从而获得第一状态和第二状态下分别与点μ₁或μ₂的欧式空间距离最近的点，称作第一邻居点，空间距离第二近的点称作第二邻居点，以此类推，空间距离第W²近的点称作第W²邻居点，将第一状态下第m次重复实验中的第k邻居点指定为P_mk ⁽¹⁾、将第二状态下第m次重复实验中的第k邻居点被记做P_mk ⁽²⁾，其中k<＝W²；

v)基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k，从而获得W²个判定标准p_k；和/或

基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中μ₂到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k’，从而获得W²个判定标准p_m’；

vi)利用rOP方法统计上述(W)²个p_m值和/或(W)²个p_m’值，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

21、根据实施方案20所述的系统，其中所述的空间泊松分布模型的分析如下：

将点P_mk ⁽¹⁾与P_mk ⁽²⁾的接触强度分别使用

和

表示，其中

如果比值

明显偏离整体分布，可认为在μ₁处第k邻居点的强度发生变化。在零假设

的情况下，比率R_k(μ₁)分别遵循具有2n₁k和2n₂k自由度的Fisher分布。

22、根据实施方案20或21所述的系统，其中所述rOP检验基于以下假设：

其中，θ_k是检验k的有效大小.S_r是p值的第r阶统计量，并且

S_r～Beta(r，2W²-r+1)，其中r值为小于p*W²的最大整数值，其中p∈(0,1]，具体含义表示两个条件之间一个染色质相互作用被视为差异性染色质相互作用时，所含有的产生显著变化的第k邻居的百分比。

需要明确的是，此处仅仅举例说明了要求保护的一些具体实施方案，其中一个或更多个技术方案中所记载的技术特征可以与任意的一个或更多个技术方案相组合，这些经组合而得到的技术方案也在本申请保护范围内，就像这些经组合而得到的技术方案已经在本发明公开内容中具体记载一样。

附图说明

图1显示了模型的过程，其中图A表明一段DNA在三维空间中的构象可能通过Hi-C的接触频率矩阵来表示。图B：统计人K细胞22号染色质体中，相隔不同距离(bin)的位置上染色质接触频率差值的方差。分别在5k分辨率和10k分辨率下(因此染色质区段的大小分别为5kb和10kb)进行了分析。图C：将染色质相互作用矩阵的位置(i，j)和该位置上的染色质相互作用强度转化为一个三维空间坐标系。图D：利用空间泊松过程计算第K邻居点的分布。

图2显示了算法参数的性能分析。图A：在模拟染色质相互作用差异的数据中，测试使用不同最近邻(例如k＝1表示使用排名第一的最近邻，k＝3表示使用排名第三的最近邻)在分析染色质相互作用差异中的效果。Prob signif表示使用所有最近邻后分析得到的结果。图B：在模拟染色质相互作用差异的数据中，比较使用不同最近邻进行泊松分布估计后进行Fisher检验的P-values的分布。图C：比较使用不同最近邻获得的ROC曲线。

图3显示了在模拟染色质相互作用差异的数据中，比较在不同窗口宽度W和不同r-th ordered p-value的条件下，本发明方案与目前主流方法edgeR在分析不同变化差异倍数的染色质相互作用时的AUC值。

图4显示了在比较GM12878和K562细胞系的染色质相互作用数据中，通过基因组特征评价本方案与目前主要方案edgeR的性能差别。A：本方案的能够发现超过1万6千个以上的染色质相互作用差异。edgeR能够发现约1万2千染色质相互作用差异。B：基因组研究发现，大部分的染色质相互作用的距离是在1M bp范围内容，超过1M bp以上的信号很可能是噪声。本方案相对edgeR所发现的染色质相互作用最大部分都落在1M bp范围内。C：基因组研究发现，大部分的染色质相互作用都是在染色质拓扑相关结构域(TAD)发生。本方案相对edgeR所发现的染色质相互作用差异最大部分都落在染色质拓扑相关结构域范围内。D：基因组研究发现，大部分的染色质相互作用与CTCF在基因组上的结合位点相关。本方案相对edgeR所发现的染色质相互作用差异更靠近CTCF在基因组上的结合位点。

图5显示了在比较GM12878和K562细胞系的染色质相互作用数据中，通过转录组特征评价本方案与目前主要方案edgeR的性能差别。A：染色质相互作用差异往往与基因表达密切相关。使用本方法发现的染色质相互作用差异的基因组位置上，相对于edgeR方法出现了更大差异倍数的基因表达变化。B，C，D：染色质相互作用差异往往与组蛋白H3K3me3修饰的位置、转录因子Pol II在基因组上的结合位置、EP300在基因组上的结合位置密切相关。本方案相对edgeR所发现的染色质相互作用差异更靠近上述信号在基因组上的位置。

图6显示了基因alpha-globin位置分析使用本方法分析在K562细胞(左图)和GM12878细胞(右图)中染色质相互作用的差异。通过本方法，我们可以容易的发现在K562细胞中出现了染色质相互作用，而这种染色质相互作用在GM12878细胞中消失了。

具体实施方式

还可进一步通过实施例来理解本发明，然而，要理解的是，这些实施例不限制本发明。现在已知的或进一步开发的本发明的变化被认为落入本文中描述的和以下要求保护的本发明范围之内。

样品

本发明所使用的样品可以是任何可被分析的对象，只要该分析的对象内部包含染色质以及基因的表达产物(例如mRNA和/或蛋白质)，样品可以是真核细胞，例如动物细胞、植物细胞、真菌细胞等，有时候也可以包括细胞的裂解物。

通过某种特定的外加诱导或者内部的自然过程而导致样品的性质或形态等发生改变，这个过程被称作状态转变。例如化学试剂的诱导分化、物理刺激或细胞在自然生理过程中的分化，如响应外部激素或其他信号分子或者细胞内部基因或蛋白的作用，而引发的细胞自然分化过程。

术语“处于第一状态的样品”和“处于第二状态的样品”是指经过状态转变后得到的两种不同状态的样品，或者是两种具有不同性质的独立样品，其中在一些实施方式中，“处于第一状态的样品”是状态转变前的样品，而“处于第二状态的样品”是经过状态转变后的样品。

染色质相互作用分析模型

理论基础

Hi-C是检测染色质空间构象的关键技术，基于Hi-C技术又演变出了多种染色质相互作用的分析技术，例如in situ Hi-C或BL-Hi-C等，通过Hi-C可以产生全基因组范围存在的大规模染色质相互作用的数据，现有技术中所有可用于染色质空间构象、染色质相互作用分析的方法均可用于本发明的方法中，以产生染色质相互作用的数据。

在本发明中，染色质相互作用的基本分析单元叫做染色质区段，在本发明中也被叫做bin，被染色质相互作用分析(如Hi-C)发现具有相互接触信号的两个(有时候是多个)染色质区段被称为染色质区段对，如“作用染色质区段对”或“邻居染色质区段对”，染色质区段对的大小主要由染色质相互作用分析方法的分辨率确定，例如，如果Hi-C的分辨率为25kb，那么区段大小可以选择25kb；如果Hi-C分辨率是5kb，染色质区段的大小就可以选择5kb，当然区段也可以选择分辨率的整倍数。

具体的，本发明的染色质区段的大小可以在0.1-100kb之间，例如0.1kb、0.2kb、0.3kb、0.4kb、0.5kb、0.6kb、0.7kb、0.8kb、0.9kb、1kb、2kb、3kb、4kb、5kb、6kb、7kb、8kb、9kb、10kb、20kb、30kb、40kb、50kb、60kb、70kb、80kb、90kb或100kb。

术语“染色质接触”是指部分染色质区域或区段之间由于空间折叠相互相靠近，从而在染色质相互作用的分析如Hi-C中，显示出相互相互作用的信号。这些信号经过分析转换为频率值后，被称为“接触强度”、“染色质接触强度”、“接触频率”或者“染色质接触频率”。其中，在对染色质相互作用的数据进行分析时，当某些相互靠近的染色质区段的接触频率超出了基准的水平或者根据本领域的普通技术知识建立的标准，即可被认为所述染色指区段之间产生了“染色质相互作用”。

术语“染色质相互作用差异”是指在不同状态的样品中，全基因组范围内或者感兴趣的区域内染色质相互作用发生改变，也可以被称为“差异性染色质相互作用”，或者“染色质相互作用变化”，在本发明中，也被简称为DCI，有时候也被简称为DI。

现有技术中基于Hi-C数据来检测染色质相互作用差异时，在建立算法的统计模型时，通常会假设特定位置发生的染色质相互作用是孤立的事件，并不会受到周边序列的影响。具体而言，假设染色质DNA位置i和位置j出现相互作用，则i位置的邻居位置(例如i+1位置，i-1位置)与j位置的邻居位置(例如j+1位置，j-1位置)的相互作用频率是不受到i位置和j位置发生相互作用影响的，同时i位置和j位置发生相互作用的频率，也不会受到i位置的邻居位置与j位置的邻居位置之间相互作用的影响。因此一个染色质位置(或DNA位置、基因组位置)与其他位置，哪怕是临近位置也是孤立存在的。

申请人经过长期的研究认为，由于染色质包含了连续存在的基因组，一个位置与邻近位置在染色质相互作用频率方面显然存在着必然的关联，在真实情况中，申请人认为如果染色质位置i和位置j出现相互作用，i位置的邻居位置(例如i+1位置，i-1位置)与j位置的邻居位置的相互作用概率也会随着i位置和j位置发生相互作用频率(概率)而增大。因此现有技术的假设本身存在很大的漏洞，从而给准确分析染色质相互作用带来很大的障碍。

模型

将完整的染色体(或者叫做染色质，在本发明中具有相同的含义)或者染色质的部分区域按一定的大小分割成若干个染色质区段(或者叫做染色体区段，bin)，按照从上游到下游的排列，依次将其命名为位置1、2、3……i-1、i、i+1、……j-1、j、j+1……，其中位置1就是表示从染色质起始位置开始的第一个染色质区段，位置2表示第2个染色质区段……，以此类推，位置i即表示第i个染色质区段，位置j表示第j个染色质区段，当染色质的位置i与位置j在空间上临近发生相互作用时，即可以使用坐标(i,j)表示一个染色质相互作用的位点。

其中上述染色质区段(bin)的大小与染色质相互作用分析方法的分辨率有关，通常而言，高分辨率获得的染色质相互作用更加精细(图1B)，在高分辨率下，可以使用更小的bin，如1-10kb，具体如1kb、2kb、3kb、4kb、5kb、6kb、7kb、8kb、9kb或10kb。但在更低的分辨率下，本发明的方法仍然使用，例如bin可以取10kb-100kb的范围内。

现有技术常见的染色质相互作用差异的分析方法通常只考虑发生相互作用的两个染色质区段本身，但是发明人发现该分析会导致很大程度的遗漏，真实的相互作用信息需要结合发生相互作用的染色指区段本身以及其周边邻居之间的接触频率信息才能够判定。因此发明人进一步将基因组的线性信息转换为一个二维平面窗口，i和j分别表示横、纵坐标值，邻居点的横坐标包括(i+1)、(i+2)……(i+w)……，(i-1)、(i-2)……(i-w)……，纵坐标包括(j+1)、(j+2)……(j+w)……，(j-1)、(j-2)……(j-w)……，这些横纵坐标的交叉处就表示了窗口内所有染色质相互作用的发生位置，设定该平面窗口的宽度为W，则窗口内一共包括W²个相互作用的位置。

窗口的宽度W通常会根据计算量和实际的需求来决定，宽度W并非越大越好，因为如果过大的W值会带来太多的不相关噪声，影响真实相互作用的判断。在本发明的一个实施方式中，W在3-19之间的奇数，W为3-19之间的奇数，例如3、5、7、9、11、13、15、17、19，优选3、5或7，在一个实施方式中，W为3，在另一个实施方式中W为5。需要说明的是，W选择奇数是为了分析和计算的方便，即选取相互作用位置两边对称的区段，同时也可以选择偶数，如2、4、6、8、10，此时只需要在后续分析和计算中进行调整即可。

二维坐标仅能反映染色质相互作用的发生位置，并不能体现出相互作用的强度或频率，因此可以进一步将上述W²个相互作用的位置及其的相互作用频率值结合，形成一个染色质相互作用频率矩阵窗口，就可以表现出更多的信息，图1A即展示了在一个具体的实施方式中，基于特定位置(i,j)及周边位置以及相应的频率值所建立的染色质相互作用频率矩阵窗口中，从中可以看出，在任意一个给定的具体位置如(i,j)，就可以获得染色质i位置和j位置的相互作用强度或频率，用Z(i，j)表示，从灰度即可表明相互作用的强弱。

将[i，j，Z(i，j)]表示为三维空间中的一个特定点，经过n_c次重复试验后，可以得到n_c个(i,j)处的相互作用频率值，进一步可以计算在第一种状态下(i，j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)n⁽¹⁾)，并在前述三维坐标系中指定一个新的点值μ₁,其坐标为[i，j，mean(Z(i，j)n⁽¹⁾)]，μ₁是虚拟出来的点，通过将窗口中W²个位置与μ₁点之间计算欧式空间距离，从而获得距离最近的第一邻居点、距离第二近的第二邻居点…距离第k近的第k邻居点，其中k<＝W²。

设定第一种状态下第m次重复试验(m<＝n_c)中的第k邻居点为p_mk ⁽¹⁾，由于实验重复n_c次，因此可获得n_c个第k邻居点，这些邻居点符合空间泊松分布模型(可以想象为一个由p_mk ⁽¹⁾的最大值作为半径的球体，这些第k邻居点在球体内部)，基于上述模型可以获得在第一种状态下观察到第k邻居的概率值；类似的，基于第二种状态下每一次重复试验中μ₁到第k邻居点的欧氏空间距离，获得在第二种状态下观察到第k邻居的概率值，利用Fisher检验比较上述概率值获得判定标准p_mk，一共可以获得W²个p_mk值；类似的，如果以μ₂为参照点，同样也可获得W²个p_mk值，为了方便区分，基于μ₂获得的p值被称为p_mk’，因此本发明的方法最多可以获得2*W²个p值(p_mk+p_mk’)。

术语“欧式空间距离”，也叫做“欧式距离”、“欧几里得度量(euclidean metric)”，其是一个通常采用的距离定义，指在m维空间中两个点之间的真实距离，或者向量的自然长度(即该点到原点的距离)。在三维空间中的欧氏距离就是两点之间的实际距离。在本发明中，计算两组染色质区段对的欧式距离时，其位置(i,j)加上染色质接触的频率Z(i,j)，实际上就形成了两组三维坐标，从而计算其欧式空间距离或欧式距离。

本发明的一个实施方案中，后续计算是基于W²个p_mk值；在一个实施方案中，后续计算是基于W²个p_mk’值；在又一个实施方案中，后续计算是基于W²个p_mk值和W²个p_mk’值。

随后使用r阶p值(rOP)统计来估计来自于2*W²检验中的r是否具有显著性。我们定义

p∈(0，1]，r表示两个条件之间的某一个染色质相互作用被视为差异性染色质相互作用时，所含有的产生显著变化的第k邻居的百分比，给定的2*W²个p值通过2*W²个检验进行估计，其中rOP统计基于以下假设：

其中，θ_k是检验k的有效大小.如果S_r是p值的第r阶统计量，则：

S_r～Beta(r，2W²-r+1)。

材料与方法

Hi-C数据

(1)Lieberman-Aiden小组发布的Hi-C数据(Rao et al.2014.A3D Map of theHuman Genome at Kilobase Resolution Reveals Principles of ChromatinLooping.Cell 159:1665–1680.)。

(2)NCBI数据库中GM12878细胞系(人B淋巴细胞)的公共数据GSM1551574和GSM1551575；以及K562细胞系(人慢性髓系白血病细胞)的公共数据GSM1551620和GSM1551623。使用juicebox工具中的VC-squared法对Hi-C矩阵进行归一化(Durand etal.,2016.Juicebox Provides a Visualization System for Hi-C Contact Maps withUnlimited Zoom.Cell Syst 3:99–101.)。对于每个不同状态的样品，使用MA-plot法进行样品间的标准化(具体描述见后文)。

基因表达数据

GSE78625(针对GM12878)和GSE78553(针对K562细胞)均来自ENCODE项目的基因表达数据。

基因集富集

使用GAGE法进行基因集的富集(Luo W et al.,2009.GAGE:generallyapplicable gene set enrichment for pathway analysis.BMC Bioinformatics 10:161.)

半方差函数计算

使用R中的gstat软件包(Pebesma 2004)进行半方差函数计算。变差函数(或半变函数)度量广泛用于地理统计中，以描述给定空间过程的空间相关程度。如果y(x)是在二维空间中定义的空间过程，所有被距离h隔开的点对之间的变化性则可以如下计算：

N(h)表示数据集中所有被距离h分隔的点，|N(h)|表示了数据集的大小。

空间过程的原则(Principles of a spatial process)

空间点过程(Spatial Point Process)是描述在d维空间中点分布的模型。如果R是一个区域，C是R的子区域，那么观察C点的概率是P(X(C)＝1|R)＝a(C)/a(R)，其中，as a(C)是C的面积。在大的参考区域R_m中，其中给定n个点，则在空间中观察到一个点的概率是：

随后

P(X(C)＝k|λ)～Poission(λ). (3)

上式也通常用于表示空间齐次泊松分布过程。

第k邻居点强度变化的显著性检验：

第k邻居点强度变化检验是基于由Burguet等人开发的k近邻密度估计器(BurguetJ,Andrey P.2014.Statistical comparison of spatial point patterns inbiological imaging.Stat Comp Spat point patterns Biol imaging；Burguet J etal.2009.Three‐dimensional statistical modeling of neuronal populations:Illustration with spatial localization of supernumerary neurons in the locuscoeruleus of quaking mutant mice.J Comp Neurol 513 VN-r，以及Burguet J,MaurinY,Andrey P.2011.A method for modeling and visualizing the three-dimensionalorganization of neuron populations from replicated data:properties,implementation and illustration.Pattern Recognit Lett.)。其主要的模型如下：

将三维空间中的Hi-C相互作用，每个相互作用可以用笛卡尔坐标(i，j，f_ij)表示，其中i和j是基因组坐标，f_ij是相互作用频率。对于给定的配对相互作用(i，j)和相互作用频率μ，在第n次重复试验中，观察到距离(i，j，μ)距离为x_ik的第k邻居点的概率是：

上式中唯一需要计算的参数是密度λ(u)(为了清楚的原因写作λ)。对于状态c进行n_c次重复的Hi-C矩阵，μ周围第k邻居点的密度，可以从上式的最大似然值估算出：

如果进行变量改变

并且替换其在等式(4)中的值，则可以得到：

替换等式(5)中的y_n，k值，则得到：

另等式7中的

可以得到:

因此，另

和

分别表示第一个状态和第二个状态下μ的第k邻居点强度，那么下式(9)中

和

的比率遵循Fisher-Snedor分布。

那么在零假设

下，在两个生物状态之间μ周围第k邻居点密度之间的比率为：

因此，给定方程10并给定置信水平α，μ周围第k邻居点密度变化的双侧p值被估计为：

整体算法流程

给定两个不同状态c₁和c₂下进行的两组Hi-C实验结果，分别进行n₁和n₂次重复。对于给定的两个基因组位置i和j之间的相互作用，令μ₁代表c₁状态下(i，j)位置的染色质相互作用频率，μ₂代表c₂状态下(i，j)位置的染色质相互作用频率，设定W代表围绕(i,j)点的周围窗口大小，W ²将包括在笛卡尔坐标系中的所有(i，j)位置附近的成对相互作用。

表示小于

的最大整数值。

使用第一状态下(i，j)的平均相互作用频率μ₁作为参考点，利用方程5、10和11，可以为所限定的窗口中的每一个相互作用关联一个p值，用于指示在两个状态下μ₁周围是否发生显著变化，从而将得到W²个p值。然后，使用第二状态下(i，j)的平均相互作用频率μ₂作为参考点，也可以计算W²个p值，用于表示μ₂周围限定窗口中相互作用的变化。因此相加将一共得到2*W ²个p值。

使用r阶p值(rOP)统计来估计来自于2*W²检验中的r是否具有显著性。我们定义

p∈(0，1]，r表示两个条件之间一个染色质相互作用被视为差异性染色质相互作用时，所含有的产生显著变化的第k邻居的百分比，给定的2*W²个p值通过2*W²个检验进行估计，其中rOP统计基于以下假设：

其中，θ_k是检验k的有效大小.如果S_r是p值的第r阶统计量，则：

S_r～Beta(r，2W²-r+1) (14)

模拟分析

以K562Hi-C强度分布图作为参考，使用负二项来模拟不同重复实验组的染色质相互靠近(作用)频率。然后随机选择10％的相互作用作为差异性的相互作用。所获得的新均值用于第二种状态下不同重复试验中差异性染色质相互作用区域的频率计数采样。

ROC计算

我们使用ROCR软件包来估计每个edgeR和本发明方法的预测性能。真是的阳性信号模拟分析中选作DI，同时算法也报告为DI的信号；假阳性信号是在模拟分析中不是DI，但由算法报告为DI的信号；假阴性信号是模拟分析中作为DI，但是由算法报告为不是DI的信号；真实的阴性信号是值模拟分析中不作为DI，同时也被算法报告为非DI的信号。

MA归一化

为了确保同一状态下的样品不同次重复之间结果的一致性，进行了样本间的归一化处理，处理过程与用于基因表达的MA-plot归一化相似。简言之，在两次重复中分别的相互作用频率为

和

的相互作用位置(i，j)，我们分别对强度(A)和比率(M)取对数如下：

然后通过在MA图中拟合loess curve来估计预期偏差M_bias。经过校正的M_correct＝M-M_bias，

和

的重新缩放值(rescaled value)的计算方法如下：

实施例1高分辨率Hi-C中的局部结构变化

由于染色质在本质上属于一种连续的聚合物，在染色质环的形成过程中，聚合物链会发生弯曲，从而使两个空间上邻近的位置i和j之间发生染色质的相互作用。当i与j之间的距离发生变化时，相邻位点(例如i-1、i+1、j-1和/或j+1)之间的距离也会发生改变。在Hi-C接触图中，这种现象表现为以成对相互作用(i，j)为中心的二维窗口内的相互作用频率的变化(图1A)。在这个窗口中存在的所有其它相互作用被称为(i，j)的“邻居相互作用”。为了定量测量Hi-C接触图中相邻位点之间空间依赖性的程度，我们在不同的Hi-C分辨率下计算了方向变异函数(图1B)。变异函数广泛用于空间统计中，用来描述给定空间过程的空间相关程度。在用于分析Hi-C接触图时，染色质相互作用的关系即通过变差函数显示，用于测量所有位置(i，j)和与之在水平、垂直或对角线方向上距离为h的位置之间的平均染色质相互作用频率差异。

图1B是基于K562细胞系第22号染色体，分别以5kb和10kb分辨率(即一个染色质区段或者bin是5kb或者10kb)归一化的Hi-C接触图(Rao等2014)中，距离为30个染色质区段时得到的方向变异图。从中可以看出，在图中所有方向上，短距离分隔的相邻位置之间具有较小的变化，同时，接触频率沿着对角线的变化比其他方向要小得多。这表明相互作用(i，j)的“影响区域”限于其临近位置之间，并且与远程的相互作用关联不大。

以上结果表明，在检测染色质相互作用差异时，需要考虑到相邻位点之间的空间相关性。

实施例2利用空间泊松过程评估染色质相互作用差异

将Hi-C作用图的染色质相互作用频率转换至三维参数空间中，其中x轴和y轴代表基因组位置坐标，z轴代表作用频率，则相互作用强度较大的区域将形成山型结构(图1C)。当存在差异相互作用时，基于本发明提出的空间依赖性假设，“山”的形态将会出现显著的变化(图1C)。为了估算两种状态下“山”的形态变化，将第一个状态下“山”尖的三维位置作为参考点(高度μ即山顶端的相互作用频率)，然后计算第一个和第二个状态下μ附近的邻居相互作用频率的变化(图1D)。

基于两个不同状态(条件)的细胞获得了Hi-C强度图，具体而言，这两个不同状态的细胞分别为：GM12878细胞系(人B淋巴细胞)以及K562细胞系(人慢性髓系白血病细胞)，Hi-C数据来自与公共数据，具体为来自NCBI数据库中GM12878细胞系(人B淋巴细胞)的公共数据GSM1551574和GSM1551575；以及K562细胞系(人慢性髓系白血病细胞)的公共数据GSM1551620和GSM1551623。每个状态都进行n_c次复制(n_c∈{1,2})。设(i，j)是我们所关注的是否发生改变的染色质相互作用，W是以它为中心的窗口(例如正方形窗口)宽度。例如，当窗口宽度为3时，将包括以位置{i-1，i，i+1}与{j-1，j，j+1}之间的笛卡尔乘积结果作为坐标的共九对相互作用。还定义μ₁为条件1下n₁次重复实验中位置(i，j)的平均相互作用频率，μ₂为条件2下n₂次重复实验中位置(i，j)的平均相互作用频率。

在这个模型中，假设参考频率为μ₁(或μ₂)的相互作用(i，j)的邻居接触频率(高度)的分布遵循齐次空间泊松分布过程。在非结构性变化的情况下，可预计可从两个生物条件下相同的空间过程中采集到染色质相互作用(i，j)的邻居相互作用频率。基于此零假设，预期在与参考点(i，j，μ₁)的距离为x的位置能够观察到第k邻居(KNN)的概率在两个生物条件之间是相似的。结果显示这个概率只取决于参考点(i，j，μ₁)周围的KNN的密度(见方法部分)。因此，在大小为W的窗口中，如果大部分KNN显示出了密度的显著变化，那么就可以认为位置(i，j)具有显著性差异的染色质相互作用。

以(i，j，μ₁)为参考点(图1D)，对每个条件的每次重复实验，可以根据与(i，j，μ₁)点的距离对(i，j)周围的窗口中W²个染色质相互作用进行排序。令

表示在条件c的第n次重复中，与参考点距离第k近的邻居。对于一个固定的k(例如第一邻居)，可以使用由Burguet开发的点强度估计器来估计第一个条件下μ₁周围KNNs的强度

以及第二个条件下μ₁周围KNNs的强度

这些强度分别使用

和

表示(分析细节参见方法部分)。如果比值

明显偏离整体，可认为在μ₁处KNN的强度发生变化。在零假设

的情况下，比率R_k(μ₁)分别遵循具有2n₁k和2n₂k自由度的Fisher分布(见“材料与方法”)。

使用(i，j，μ₁)以及(i，j，μ₂)作为参考，对每个点计算2*W²个p值。然后将这些p值使用r阶p值统计(rOP)(Song和Tseng，2014)合并，并使用Benjamini-Hochberg程序(Benjamini和Hochberg，1995)进行校正(参见“材料与方法”)。

实施例3相邻的相互作用对Hi-C作用频率变化具有不同的敏感度

为了检测上述模型的性能，同时检查每个k近邻对最终决策的贡献，我们模拟了两个Hi-C条件，每个条件都有两个重复，并使用高斯核函数来平滑DCIs区域周围的频率。随后为每一个邻居赋予相同的权重。设sig是一个2*W²的向量，如果p值不显著，则sig(i)＝0，否则等于1。如果至少有一半的最近邻居在密度μ₁和μ₂周围的密度发生显著变化，则位置(i，j)具有差异的染色质相互作用：

为了评估每个KNN的效果，为每个在两条件之间第k邻居相互作用频率具有显著差异的KNN进行了制图(k∈[1,9])。结果显示，最远的邻居对变异性更加敏感，这些邻居报告了很多显著性强度变化，相对而言，最近的邻居则对变异性较不敏感。

对每个KNN获得的p值分布进行分析(图2B，其中1-9使用μ₁作为参考，10-18使用μ₁作为参考)也验证了类似结论：远距离邻居有许多p值位于零假设接受区

以外的末端，而最近的邻居有更多的p值在位于零假设接受区域内。事实上，这种灵敏度反映了Fisher分布累积密度函数(CDF)的收敛速度，如图2C所示，自由度越小，CDF收敛到1的速度就越慢。因此，对于较小的k值，如果使比率Rk(u)出现显著性，它需要非常小

或非常大

相比之下，远距离的邻居只需要在零假设拒绝区域中具有小的可变性。

可以采用不同的策略来组合每个成对相互作用(i，j)所得到的2*W²个p值。如Fisher's组合概率检验(Fisher R.1925.Statistical methods for research workers.)和Stouffer's Z检验(Riley JW et al.1949.The American Soldier:Adjustment DuringArmy Life.Am Sociol Rev 14:557.)或其加权变体，然而,上述方法中如果2*W ²测试中有一个是非空的，则会得出显著性的p值，这使得它们对最远的KNN变化极其敏感。诸如maxP(Wilkinson B.1951.A statistical consideration in psychologicalresearch.Psychol Bull 48:156–8.)等其它方法则过于严格，因为其只考虑2*W ²测试中最高阶的p值，因而会漏掉很多差异性相互作用。理想情况下，应该采用一个加权的p值替代方案，然而这往往对计算的要求很高。一个很好的折衷是使用r阶排序的p值(rOP)方法。如果至少一组检验的r检验是显著的，则rOP认为其具有显著性差异。在零假设下，使用自由度为α＝r和β＝2W²-r+1的β分布，对排序的2*W²个p值中的r阶统计量进行rOP检验。在本模型中，r值为小于p*W²的最大整数值，(p∈(0,1])，表示两个条件之间一个染色质相互作用被视为DCI时，所具有的显著变化的KNNs的百分比。

实施例4本发明方法与edgeR基于模拟数据的比较

现有技术中的edgeR方法即假定染色质片段之间的相互作用是独立存在的，被作为代表与本发明的方法(如实施例1-3中所披露的完整流程)进行比较。模拟进行了在两个条件下的Hi-C实验，各自进行两次重复。使用5kb分辨率的K562Hi-C的强度信号，其中，其negative binominal均值等于K562矩阵中的成对相互作用时，被认为是无差异相互作用并富集出来；而其均值等于其对应的K562Hi-C接触图中的成对相互作用的倍数变化时，被认为是差异性相互作用。通过尽可能通过选择尽可能少量的相互作用作为DCI(大约1％)，从而尽可能使模拟的DCIs更加稀疏。在这些DCI中，在一个给定的倍数变化下，大约40％具有相互作用技术的增加，在差异作用染色质区段周围应用高斯平滑器来模拟邻居变化的影响。

基于上述设置，使用曲线下面积(AUC)值来比较在给定不同窗口大小(W)、显著变化的KNN(p)百分比和倍数变化(FD)值时，本发明方法相对于EdgeR的性能差别。

对于每个固定设置，模拟重复10次，结果显示在图3中，对于小于或等于5个bin(25kb)的窗口，本发明方法明显优于edgeR。不过选择差异性KNN的比例，会在一定程度上影响到本发明方法的性能。当要求所有邻居染色质区段的接触频率差异均需要达到显著变化时(即p＝1)，虽然本发明方法仍然能够有效的检测到DCI，但是所检测到的DCI比例明显下降。因此，这种严格的设置对于分析目的在于仅检测最明显的DCI时更加有效，因为有可能会遗漏一些差异并不非常显著但却比较重要的DCI；相反，如果只要求少量邻居染色质区段的相互接触达到显著性差异(例如p<0.5)，即可实现很高的检测灵敏度和稳定性。此外，当窗口的尺寸非常大(例如超过35kb时)，本发明方法的性能出现一定程度的下降，其可能的原因在于如果窗口尺寸过大，相当于在更大范围内分析染色质相互作用改变，相当于引入了过多的无关信息，掩盖了局部染色质相互作用的信息，事实上导致染色质相互作用分析的分辨率下降。

此外，虽然edgeR在高倍变化区域中的效果类似于本发明的方法，例如，在允许DCI区段具有10或更多的倍数变化的情况下，edgeR和本发明方法的效果类似，但是在低倍变化区域的结果表明，edgeR遗漏了很多位点本身显示较低的染色质相互作用强度的变化、但在其附近具有显著结构变化的DCI，而这些DCI可以被本发明的方法有效检测到(图3)。

实施例5本发明方法与edgeR基于实际Hi-C数据的比较

将K562和GM12878细胞系的5kb分辨率Hi-C接触图(Rao et al.2014.A 3D Map ofthe Human Genome at Kilobase Resolution Reveals Principles of ChromatinLooping.Cell 159:1665–1680.)进行比较，对于每个细胞系进行考两次重复，窗口宽度W选择3或5，数据用knight-Ruiz矩阵平衡算法进行归一化处理(参考文献同上)，此外还使用MA-plot法(见“材料与方法”部分)对每个细胞系数据的不同的重复试验进行归一化。

本发明方法所检测到的染色质相互作用差异是edgeR法所检测的1.6倍(图4A)。

已知大多数染色质结构变化倾向于发生在拓扑相关域(TAD)内，因此接下来分析了两种方法所获得的DCI的跨度，结果发现本发明方法所检测到的DCI具有小于1Mb(平均跨度58229.4kb)，而edgeR法检测到的大约20％的DCI跨度大于1Mb，平均跨度为107,555kb(图4B)。进一步计算了DCI位于TAD内部的比例(图4C)，结果显示本发明方法检测到的DCI中约70％位于TAD内，而对于edgeR，只有约20％的DCI位于TAD内。另外，由于CTCF对形成和维持染色质整体结构具有关键作用，初步预期DCI应位于差异性CTCF峰值的附近。于是接下来计算了edgeR和本发明方法检测到的DCI到差异性CTCF峰值的距离。与edgeR相比，本发明方法检测的DCI有更大比例位于CTCF峰附近的小于100bp的范围。所有上述结果都表明了本发明方法比edgeR更加可靠。

染色质拓扑结构的改变对增强子和启动子之间的关联产生很大影响，从而改变基因表达，随后将位于DCI附近的基因根据其表达倍数变化(FC)分为两类：FC<＝2(在两个细胞系之间表达量不具有显著差异的基因)；FC>2(在两个细胞系之间表达量具有显著差异的基因)。对于本发明方法，高达71.46％的DCI相关基因的表达量具有显著差异(FC>＝2)，而对于edgeR，只有约50.63％的DCI相关基因在细胞系之间出现表达差异(图5A)。同时，与前述结果一致，通过本发明方法检测到的DCI更接近具有表达活性的基因信号，如H3K4me3、Pol II结合位点和EP300(图5B-D)。以上结果表明，本发明方法相比于比edgeR能够更多的检测到与实际功能相关的DCI。

实施例6染色质结构在K562细胞分化中的作用

染色质环的作用在经典的α-珠蛋白基因座中得到了很好的表征，该基因座在红系细胞(K562)中具有排他性表达，同时在淋巴母细胞中(GM12878)的沉默表达(图6)。这个基因座上的基因受到位于α-球蛋白基因上游30～60kb的远端DNase I超敏位点(HS)的调控，并且已知所述基因的沉默是由于缺少增强子与启动子的相互作用而导致的。

如图6所示，通过设定p值cut-off为1e-4对结果进行标记，发现本发明方法所有报道的具有差异的染色质相互作用都位于与CTCF蛋白结合出现差异的位置附近。在这个区域的ChromHMM轨迹也表明此处的染色质状态变化很大，从强增强子状态(在ChromHMM轨迹中以橙色显示)转变为异染色质状态(在ChromHMM轨迹中以灰色显示)。

Claims (20)

1.一种用于分析不同状态的样品中染色质相互作用差异的方法，其包括下列步骤：

(1)分别对处于第一状态和第二状态的样品进行染色质相互作用的分析；

(2)根据作用染色质区段对在第一状态和第二状态的样品间的接触强度差异以及邻居染色质区段对之间的接触强度差异，分析得到在不同状态的样品中具有差异的染色质相互作用，

所述作用染色质区段对是指经过染色质相互作用分析被鉴定实际发生染色质相互作用的染色质区段对，其包括至少两个作用染色质区段；

所述邻居染色质区段对是由邻居染色质区段两两组合形成的，所述邻居染色质区段是指发生相互作用的两个染色质区段各自所临近的染色质区段，由于作用染色质区段对的关系，作用染色质区段附近的染色质区段之间的空间距离也被拉近而产生染色质接触的信号；

其中所述步骤(1)的染色质相互作用分析重复进行n_c次实验，其中c＝1或2，n₁表示针对第一状态的样品重复的次数，n₂表示针对第二状态的样品重复的次数，n₁和n₂相等或不相等，n_c是1-100之间的自然数；

其中所述步骤(2)中获得染色质相互作用的差异包括下列步骤：

a)计算处于第一状态的样品中，每一个或感兴趣的作用染色质区段对在n₁次重复实验中的平均染色质作用强度；和/或处于第二状态的样品中，每一个或感兴趣的作用染色质区段对在n₂次重复实验中的平均染色质作用强度；

b)以步骤a)获得的处于第一状态的样品或处于第二状态的样品的平均染色质作用强度为基准强度，并为该基准强度赋予与所述作用染色质区段相同的空间位置信息，分别计算每一次重复实验中，处于第一状态和第二状态的样品中每一个或者感兴趣的作用染色质区段对及其邻居染色质区段对相对于基准强度的排序，所述排序综合了相对空间位置和相对染色质接触强度的信息，所述排序是基于各染色质区段对相对于基准强度的欧式距离而确定的，

c)分别获得针对每个样品在不同次的实验中具有相同排序的染色质区段对的集合，建立上述集合相对于基准强度的分布，分析每一排序下的所述分布在处于第一状态的样品和第二状态的样品之间是否具有显著性差异，所述分布是利用数学模型计算获得。

2.根据权利要求1所述的方法，其中步骤c)所述数学模型是空间泊松点过程(空间泊松分布)，基于该空间泊松点过程计算得到该分布下观察到具有该排序的染色质接触强度的概率，并分析上述概率值在第一状态的样品和第二状态的样品之间是否具有显著性差异。

3.根据权利要求2所述的方法，其中分别以处于第一状态和第二状态的样品的平均染色质作用强度为基准强度，各进行一次空间泊松点过程的计算。

4.根据权利要求2所述的方法，其中分析所述分布之间是否具有显著性差异时采用Fisher检验，获得p值，随后将由各个排序下的分布获得的p值集合进行rOP统计(r-thordered p-value statistics)，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

5.根据权利要求1所述的方法，所述邻居染色质区段对由位于作用染色质区段两侧各1-10个位置的邻居染色质区段之间两两组合形成。

6.根据权利要求2所述的方法，其中步骤(2)按照下述步骤进行：

i)为一个作用染色质区段对赋予坐标信息(i,j)，其中i表示发生相互作用的两个染色质区段中，位于基因组序列上游的染色质区段位置，j表示位于基因组序列下游的染色质区段位置；以(i,j)为平面中心建立一个二维平面窗口，宽度为W，所述窗口共包含W²个染色质区段对，用w表示位于与i或j的坐标距离，则上述窗口内各染色质区段对的横坐标可表示为(i+w)、i、和(i-w)，纵坐标为(j+w)、j和(j-w)，除了中心位置(i,j)外，其他位置作为邻居位置，表示为(i[w],j[w])；其中w<＝(W-1)/2，W为3-19之间的奇数；

ii)将第m次重复试验中(i,j)位置的染色质相互作用频率记做Z(i，j)_m，其中m是自然数，并且m<＝n_c；将位置坐标(i,j)以及该位置的频率值转换为一个三维坐标系的点(i，j，Z(i，j)_m)；同时将窗口内其他位置的染色质接触强度记做Z(i[w],j[w])_m，进一步结合位置坐标将染色质区段对转换为多个三维坐标系的点(i[w]，j[w]，Z(i[w],j[w])_m)；

iii)计算在处于第一种状态的样品中(i，j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽¹⁾)，并在前述三维坐标系中指定一个新的点值μ₁,其坐标为[i，j，mean(Z(i，j)_n ⁽¹⁾)]；和/或

计算处于第二种状态下的样品中(i,j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽²⁾)，并在前述三维坐标系中指定一个新的点值μ₂，其坐标为[i，j，mean(Z(i，j)_n ⁽²⁾)]；

iv)分别计算处于第一种状态和第二次状态下的样品的第m次重复实验中，包含中心位置(i，j，Z(i，j)_m)以及邻居位置(i[w]，j[w]，Z(i[w],j[w])_m)在内的点到μ₁或μ₂的欧式空间距离，从而获得第一状态和第二状态下分别与点μ₁或μ₂的欧式空间距离最近的点，称作第一邻居点，空间距离第二近的点称作第二邻居点，以此类推，空间距离第W²近的点称作第W²邻居点，将第一状态下第m次重复实验中的第k邻居点指定为P_mk ⁽¹⁾、将第二状态下第m次重复实验中的第k邻居点被记做P_mk ⁽²⁾，其中k<＝W²；

v)基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k，从而获得W²个判定标准p_k；和/或

基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中μ₂到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k ^’，从而获得W²个判定标准p_m ^’；

vi)利用rOP统计上述(W)²个p_m值和/或(W)²个p_m ^’值，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

7.根据权利要求6所述的方法，其中所述的空间泊松分布模型的分析如下：

将点P_mk ⁽¹⁾与P_mk ⁽²⁾的接触强度分别使用

和

表示，其中

如果比值

明显偏离整体分布，可认为在μ₁处第k邻居点的强度发生变化；在零假设

的情况下，比率R_k(μ₁)分别遵循具有2n₁k和2n₂k自由度的Fisher分布。

8.根据权利要求6所述的方法，其中所述rOP检验基于以下假设：

其中，θ_k是检验k的有效大小.S_r是p值的第r阶统计量，并且

S_r～Beta(r，2W²-r+1)，其中r值为小于p*W²的最大整数值，其中p∈(0,1]，具体含义表示两个条件之间一个染色质相互作用被视为差异性染色质相互作用时，所含有的产生显著变化的第k邻居的百分比。

9.一种鉴定调控染色质相互作用的试剂的方法，其包括将使样本与一种或多种试剂接触，利用权利要求1至8中任一项所述的方法分析染色质相互作用差异，以及

鉴定相比于不添加调控试剂的对照组能够改变相互作用的试剂。

10.一种分析细胞分化或发育过程中遗传物质高级结构改变的方法，其包括权利要求1-8任一项所述的步骤。

11.一种鉴定染色质结构变异的方法，其包括权利要求1至8中任一项所述的步骤。

12.一种鉴定能够调控遗传物质高级结构或引起染色质结构变异的调控试剂的方法，其包括将使样本与一种或多种试剂接触，利用权利要求1至8中任一项所述的方法分析染色质相互作用差异，以及

鉴定相比于不添加调控试剂的对照组能够改变相互作用的试剂。

13.一种染色质相互作用差异的鉴定系统，所述系统包括染色质相互作用实验模块和数据分析模块，其中所述染色质相互作用实验模块用于分别对处于第一状态和第二状态的样品中所存在染色质相互作用进行分析；

所述数据分析模块用于根据作用染色质区段对在第一状态和第二状态的样品间的接触强度差异以及邻居染色质区段对之间的接触强度差异，分析得到在不同状态的样品中具有差异的染色质相互作用；

所述作用染色质区段对是指经过染色质相互作用分析被鉴定实际发生染色质相互作用的染色质区段对，其包括至少两个作用染色质区段；

所述邻居染色质区段对是由邻居染色质区段两两组合形成的，所述邻居染色质区段是指发生相互作用的两个染色质区段各自所临近的染色质区段，由于作用染色质区段对的关系，作用染色质区段附近的染色质区段之间的空间距离也被拉近而产生染色质接触的信号；

其中染色质相互作用的分析重复进行n_c次实验，其中c＝1或2，n₁表示针对第一状态的样品重复的次数，n₂表示针对第二状态的样品重复的次数，n₁和n₂相等或不相等，n_c是1-100之间的自然数；

其中数据分析模块获得染色质相互作用的差异包括下列步骤：

a)计算处于第一状态的样品中，每一个或感兴趣的作用染色质区段对在n₁次重复实验中的平均染色质作用强度；和/或处于第二状态的样品中，每一个或感兴趣的作用染色质区段对在n₂次重复实验中的平均染色质作用强度；

b)以步骤a)获得的处于第一状态的样品或处于第二状态的样品的平均染色质作用强度为基准强度，并为该基准强度赋予与所述作用染色质区段相同的空间位置信息，分别计算每一次重复实验中，处于第一状态和第二状态的样品中每一个或者感兴趣的作用染色质区段对及其邻居染色质区段对相对于参比强度的排序，所述排序综合了相对空间位置和相对染色质接触强度的信息，所述排序是基于各染色质区段对相对于基准强度的欧式距离而确定的，

c)分别获得针对每个样品在不同次的实验中具有相同排序的染色质区段对的集合，建立上述集合相对于基准强度的分布，分析每一排序下的所述分布在处于第一状态的样品和第二状态的样品之间是否具有显著性差异，所述分布是利用数学模型计算获得。

14.根据权利要求13所述的系统，其中步骤c)所述数学模型是空间泊松点过程(空间泊松分布)，基于该空间泊松点过程计算得到该分布下观察到具有该排序的染色质接触强度的概率，并分析上述概率值在第一状态的样品和第二状态的样品之间是否具有显著性差异。

15.根据权利要求14所述的系统，其中分别以处于第一状态和第二状态的样品的平均染色质作用强度为基准强度，各进行一次空间泊松点过程的计算。

16.根据权利要求14所述的系统，其中分析所述分布之间是否具有显著性差异时采用Fisher检验，获得p值，随后将由各个排序下的分布获得的p值集合进行rOP统计(r-thordered p-value statistics)，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

17.根据权利要求13所述的系统，所述邻居染色质区段对由位于作用染色质区段两侧各1-10个位置的邻居染色质区段之间两两组合形成。

18.根据权利要求13至17中任一项所述的系统，其中步骤(2)按照下述步骤进行：

i)为一个作用染色质区段对赋予坐标信息(i,j)，其中i表示发生相互作用的两个染色质区段中，位于基因组序列上游的染色质区段位置，j表示位于基因组序列下游的染色质区段位置；以(i,j)为平面中心建立一个二维平面窗口，宽度为W，所述窗口共包含W²个染色质区段对，用w表示位于与i或j的坐标距离，则上述窗口内各染色质区段对的横坐标可表示为(i+w)、i、和(i-w)，纵坐标为(j+w)、j和(j-w)，除了中心位置(i,j)外，其他位置作为邻居位置，表示为(i[w],j[w])；其中w<＝(W-1)/2，W为3-19之间的奇数；

ii)将第m次重复试验中(i,j)位置的染色质相互作用频率记做Z(i，j)_m，其中m是自然数，并且m<＝n_c；将位置坐标(i,j)以及该位置的频率值转换为一个三维坐标系的点(i，j，Z(i，j)_m)；同时将窗口内其他位置的染色质接触强度记做Z(i[w],j[w])_m，进一步结合位置坐标将染色质区段对转换为多个三维坐标系的点(i[w]，j[w]，Z(i[w],j[w])_m)；

iii)计算在处于第一种状态的样品中(i，j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽¹⁾)，并在前述三维坐标系中指定一个新的点值μ₁,其坐标为[i，j，mean(Z(i，j)_n ⁽¹⁾)]；和/或

计算处于第二种状态下的样品中(i,j)位置上n_c次重复试验中得到的相互作用频率的平均值，记做mean(Z(i，j)_n ⁽²⁾)，并在前述三维坐标系中指定一个新的点值μ₂，其坐标为[i，j，mean(Z(i，j)_n ⁽²⁾)]；

iv)分别计算处于第一种状态和第二次状态下的样品的第m次重复实验中，包含中心位置(i，j，Z(i，j)_m)以及邻居位置(i[w]，j[w]，Z(i[w],j[w])_m)在内的点到μ₁或μ₂的欧式空间距离，从而获得第一状态和第二状态下分别与点μ₁或μ₂的欧式空间距离最近的点，称作第一邻居点，空间距离第二近的点称作第二邻居点，以此类推，空间距离第W²近的点称作第W²邻居点，将第一状态下第m次重复实验中的第k邻居点指定为P_mk ⁽¹⁾、将第二状态下第m次重复实验中的第k邻居点被记做P_mk ⁽²⁾，其中k<＝W²；

v)基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k，从而获得W²个判定标准p_k；和/或

基于第一种状态下每一次重复试验中得到的μ₁到第k邻居点P_mk ⁽¹⁾的欧式空间距离，使用空间泊松分布模型分析获得在第一种状态下观察到第k邻居点的概率值；同时基于第二种状态下每一次重复试验中μ₂到第k邻居点P_mk ⁽²⁾的欧氏空间距离，获得在第二种状态下观察到第k邻居点的概率值，利用Fisher检验比较上述概率值获得判定标准p_k ^’，从而获得W²个判定标准p_m ^’；

vi)利用rOP统计上述(W)²个p_m值和/或(W)²个p_m ^’值，从而确定处于第一状态的样品和第二状态的样品中在每一个或者感兴趣的染色质相互作用处是否存在差异。

19.根据权利要求18所述的系统，其中所述的空间泊松分布模型的分析如下：

将点P_mk ⁽¹⁾与P_mk ⁽²⁾的接触强度分别使用

和

表示，其中

如果比值

明显偏离整体分布，可认为在μ₁处第k邻居点的强度发生变化；在零假设

的情况下，比率R_k(μ₁)分别遵循具有2n₁k和2n₂k自由度的Fisher分布。

20.根据权利要求18所述的系统，其中所述rOP检验基于以下假设：

其中，θ_k是检验k的有效大小.S_r是p值的第r阶统计量，并且

S_r～Beta(r，2W²-r+1)，其中r值为小于p*W²的最大整数值，其中p∈(0,1]，具体含义表示两个条件之间一个染色质相互作用被视为差异性染色质相互作用时，所含有的产生显著变化的第k邻居的百分比。

Images