一种基于染色质随机折叠过程预测染色质域结构TAD的方法

一种基于染色质随机折叠过程预测染色质域结构tad的方法
技术领域
1.本发明属于染色质三维结构技术领域，涉及一种基于染色质随机折叠过程预测染色质域结构tad的方法。


背景技术：

2.哺乳动物的染色质线性长度可达2米，它折叠在直径为10-20微米的细胞核中形成复杂的三维(3d)结构。该3d结构决定了特定染色质片段之间的空间相互作用，从而启动特定的生命进程。域结构(topologically associating domains,tads)是染色质3d结构的组成单元，它代表一段染色质区域，该区域中的染色质片段倾向于只与区域内的染色质片段在空间上相互接触，与区域外的染色质片段绝缘。因此，tads对于理解基因调控、细胞分化、细胞病变等生物学过程至关重要。
3.实验上通常采用高通量染色体构象捕获技术(high-throughput chromosome conformation capture,hi-c)来获取tads的信息。该技术将多个细胞中的染色质同时冻结、交联、剪切，然后进行测序，得到细胞中染色质片段之间空间接触的平均信息。该方法的优点是，可以在全基因尺度上获取tads信息。缺点是，成本高、难以广泛应用。
4.由于染色质本质上是柔性高分子链，计算模拟方面可以通过分子动力学模拟(molecular dynamics,md)染色质3d结构来预测tads。md已用于预测染色质tads结构。现有方法通常假设染色质是一种dna密度分布均匀的高分子链。然后，在该均匀的高分子链上根据基因组测序数据假设各种定向作用，依据这些定向作用构建相应的势能函数，在此基础上建立高分子链3d结构的总势能作为md模拟的力场，进行md计算以模拟染色质3d结构预测tads。预测结果的正确性通过hi-c实验数据来验证。
5.基因组测序实验提供多种表观遗传信息，包括dna甲基化、染色质可及性、组蛋白修饰以及转录活性等。由于基因组测序实验成本低廉，已经用于临床疾病检测。因此，在基因组测序数据基础上，模拟染色质3d结构预测tads具有极大的应用前景。然而，现有的基于定向作用的预测方法通常计算量大，只能模拟包含几个基因的一段染色质区域(小于10mb)的3d结构，无法模拟全基因组(约8000mb)的3d结构，因此不能在全基因组尺度上预测tads。


技术实现要素：

6.本发明的主要目的在于提供一种基于染色质随机折叠过程预测染色质域结构tad的方法，以解决上述背景技术中提出的问题。
7.为了达到上述目的，本发明采用以下技术方案：
8.一种基于染色质随机折叠过程预测染色质域结构tad的方法，包括以下步骤：
9.(1)为待研究的染色质区间构建携带染色质可及性信息的染色质高分子链模型集合；
10.(2)模拟携带染色质可及性信息的染色质高分子链模型在微小球形空间中的随机折叠过程，计算得到染色质3d结构的集合；
11.(3)基于层级优化法得到精度高达10kb的待研究的染色质区间的3d结构集合；
12.；
13.(4)基于步骤(3)所预测染色质3d结构集合构建平均接触矩阵以显示tad结构。
14.其中，染色质视为dna密度分布不均匀的高分子链，该高分子链在细胞核中的随机折叠是染色质三维结构单元tads形成的基础。随机折叠是指高分子链在无定向作用势能作用下的折叠过程。其中无定向作用势能是指高分子链中任意两个链珠之间的相互作用势的函数形式是相同的。
15.其中，步骤(1)产生的高分子链集合中的每一个高分子链的总势能函数分为三个部分：谐振子势能e
har
、排斥势能e
wca
和空间约束势能e
sphere
；其中，谐振子势能用于保证高分子链中相邻链珠之间的连续性，排斥势能用于避免两个不相邻的链珠之间发生重叠，空间约束势用于模拟由于细胞核的狭窄所产生的空间限制；
16.本发明采用以上技术方案，提出了染色质在细胞核中的随机折叠是染色质3d结构单元tads形成的基础，构建了携带染色质可及性信息的染色质高分子链模型，用于预测染色质三维结构的实验数据(染色质可及性信息)来自基因组测序实验(atac-seq或dnase-seq实验)，有极大的实用性，可以根据现有的数据库，为各种细胞类型提供染色质3d结构中的tads信息。
17.优选的，所述步骤(1)中，构建携带染色质可及性信息的染色质高分子链模型的具体过程为：
18.步骤(11)：选取待研究的染色质区间，序列长度为l，并获取该区间对应的多细胞平均染色质可及性数据实验数据，其中，对于哺乳动物细胞，l的取值范围为：10mb《l《8000mb；
19.步骤(12)：将待研究染色质区间等分为序列长度为l的m个染色质片段，m＝l/l，并根据步骤(11)得到的染色质可及性实验数据，为每个所述染色质片段确定一个归一化的可及性信号值λ，以对每个所述染色质片段的染色质可及性性质进行定量，其中，l的取值范围为：5kb《l《5000kb；
20.步骤(13)：将每个所述染色质片段对应的粗粒化链珠的数目视为独立的随机泊松变量，其取值服从泊松分布，概率为：
[0021][0022]
其中，ki表示第i个染色质片段对应的链珠数目，λi是第i个染色质片段的染色质可及性信号值，每个链珠具有相同的直径；
[0023]
步骤(14)：由m个染色质片段构成的待研究的染色质区间的高分子链模型所包含的链珠的总数目为步骤(15)：由于第i个染色质片段的链珠数目ki是一个随机泊松变量，可以为待研究的染色质区间随机产生多个高分子链串珠模型，形成一个集合。该集合中每个模型携带确定的染色质可及性信息，它们的平均值与步骤(11)提到的多细胞平均染色质可及性实验数据相一致。
[0024]
本发明将染色质片段用粗粒化的珠子来代表，以此构建染色质高分子链模型，是一种常用方法。现有的构建染色质高分子链模型的通用方法是，用相同数目的链珠来表示
相同序列长度的染色质片段。
[0025]
本发明进一步采用以上技术特征，其优点在于，用不同数目的链珠来表示相同序列长度的染色质片段，因此可以构建携带染色质可及性信息的染色质高分子链模型。
[0026]
优选的，所述步骤(11)中，多细胞平均染色质可及性数据实验数据来自多细胞atac-seq或dnase-seq实验。
[0027]
优选的，所述步骤(1)中构建的染色质高分子链模型的精度为步骤(12)等分后染色质片段序列长度lkb。
[0028]
优选的，所述步骤(2)包括：
[0029]
步骤(21)：从不同的随机结构出发，对每条高分子链的随机折叠过程进行多次分子动力学模拟；
[0030]
步骤(22)：通过对不同高分子链和不同随机结构出发的模拟结果进行收集，得到一个待研究的染色质区间的3d结构集合。
[0031]
本发明进一步采用以上技术特征，其优点在于，模拟携带染色质可及性信息的高分子链的随机折叠过程，以此预测染色质3d结构。
[0032]
优选的，所述步骤(3)中，所述层级优化的过程如下：
[0033]
步骤(31)：构建低精度的高分子链模型，从随机结构出发进行优化；
[0034]
步骤(32)：在低精度优化结构的基础上，采用插值法产生高精度的3d结构；
[0035]
步骤(33)：以此步骤(32)得到的结构作为初始结构进行高精度的优化，得到指定的高精度结构。
[0036]
其中，优化是指进行驰豫，即以随机折叠过程运动一段时间或一定步数。
[0037]
本发明进一步采用以上技术特征，层级优化的核心思想是先构建低精度的高分子链模型，从随机结构出发进行优化。在低精度优化结构的基础上，采用插值法产生高精度的3d结构，以此结构作为初始结构进行高精度的优化，以此类推，最后得到指定的高精度结构。对于哺乳动物全基因组染色，通常初始随机结构的精度为5120-kb，优化该初始结构，以优化后的5120-kb精度的结构为基础产生2560-kb精度的初始结构。依次类推，分别产生和优化1280-kb，640-kb,320-kb,160-kb,80-kb,40-kb,20-kb以及最后10-kb精度的3d结构。
[0038]
本发明进一步采用以上技术特征，其优点在于，节省计算资源，可以快速优化精度高达10kb的全基因组哺乳动物染色质3d结构。
[0039]
优选的，所述步骤(4)包括：
[0040]
步骤(41)：如步骤(2)所述，3d结构集合中的每一个结构都是通过模拟高分子模型的随机折叠得到的。如步骤(1)所述，每个高分子模型都由m个等序列长度染色质片段组成。在任意一个3d结构中，测量3d结构中每一对染色质片段质心之间的距离，可以为3d结构构建单个距离矩阵(individual distance matrix)。
[0041]
步骤(42)：根据模型中珠链的直径设定一个距离阈值，两个片段之间的距离低于阈值时，认为该两个片段之间存在接触(contact)。据此，可以将步骤(41)中构建的单个距离矩阵转化为单个接触矩阵(individual contact matrix)。
[0042]
步骤(43)：将单个距离矩阵和单个接触矩阵在整个3d结构集合上进行平均，可得到平均距离矩阵(ensemble-average distance matrix)和平均接触矩阵(ensemble-average contact matrix)。
[0043]
步骤(44)：定义分离分数(separation score)，用以定量分析平均接触矩阵中每一个基因位点上下游染色质区域的分离程度；通过分离分数显示tads结构。
[0044]
其中，分离分数是一种用于定量分析平均接触矩阵中每一个基因位点上下游染色质区域的分离程度物理量。该分离分数数值定义为以下两组数值的差异：(1)(上游染色质区域内部相互作用的数值 下游染色质区域内部相互作用的数值)/2，(2)上游染色质区域与下游染色质区域之间相互作用的数值。其中上游染色质(或下游染色质)区间的区间长度为s，取值范围为100kb-1mb。
[0045]
本发明进一步采用以上技术特征，其优点在于，本发明预测的tads结构可以同时被单细胞结构成像(fluorescence in situ hybridization,fish)数据和多细胞hi-c构象捕获实验数据所验证。
[0046]
与现有技术相比，本发明带来了如下有益效果：
[0047]
1、本发明提出了染色质在细胞核中的随机折叠是染色质域结构tad形成的基础。
[0048]
2、本发明构建了携带染色质可及性信息的染色质高分子链模型。
[0049]
3、本发明采用层级优化法对携带染色质可及性信息的高分子链模型的随机折叠过程进行分子动力学模拟，在10kb精度下预测了哺乳动物全基因组3d结构中40％以上的tad结构。
[0050]
4、相比于以前的预测染色质tad结构的方法，本发明预测的结果首次可以同时被单细胞结构成像(fluorescence in situ hybridization,fish)数据和多细胞hi-c构象捕获实验数据验证。
[0051]
6、本发明具有极大的实用性，可以根据现有的染色质可及性数据库，为各种细胞类型预测染色质tad结构，比如健康细胞、癌细胞、病毒感染细胞等。
[0052]
7、本发明为从染色质3d结构的视角研究细胞分化过程、疾病的发生发展、以及疾病的早诊断等提供了可靠且实用的技术手段。
附图说明
[0053]
图1为实施例中从atac-seq的染色质可及性实验数据。
[0054]
图2为实施例中根据图1的数据计算得到的十种不同高分子链模型。
[0055]
图3为实施例中选取的5mb染色质区间中tad结构与多细胞hi-c实验中tad结构的对照。
[0056]
图4为实施例中全基因组染色质区间中tad结构与多细胞hi-c实验中tad结构重合的程度图。
[0057]
图5为实施例中从预测结构中计算的平均距离矩阵与单细胞fish实验数据计算得到的平均距离矩阵的对照。
[0058]
图6为实施例中图5展示的两个平均距离的相关图。
[0059]
图7为实施例中预测的细胞分化过程中关键的染色质结构重组与hi-c实验揭示的重组的对照。
具体实施方式
[0060]
为了使本技术领域的人员更好地理解本技术方案，下面将结合本技术实施例中的
附图，对本技术实施例中的技术方案进行清楚、完整地描述。
[0061]
以人类红白血病k562细胞以及从造血干细胞和祖细胞(hspc)到成熟免疫t细胞分化过程所经历的8个阶段的细胞为例，阐述基于染色质随机折叠过程预测染色质域结构tad的方法。
[0062]
对于人类红白血病k562细胞，首先选取待研究的染色质区域，它包括22条常染色体和1条x染色体。然后，从基因表达综合数据库(the gene expression omnibus，geo)获取这个染色质区域的染色质可及性数据,在本实施例中选取了atac-seq实验数据。随后，将这个包含23条染色体的染色质区域分割成连续的m(约为30万)个10-kb长度的染色质片段，并依据atac-seq实验原始数据计算每个片段的染色质可及性信号值ζ，然后对片段的可及性信号进行归一化，最后得到每个片段的染色质可及性定量数值，计算公式为：
[0063][0064]
其中，λi表示第i个染色质片段对应的归一化染色质可及性数值，ζi表示第i个染色质片段对应的染色质可及性信号数值，median(ζ)为所有染色质片段对应的染色质可及性数值的中位数。接来下，确定第i个10-kb染色质片段应该由多少个等体积的粗粒化链珠来代表。依据确定的λi值，将第i个10-kb染色质片段所对应的珠链数目ki看作一个独立的随机泊松变量，它的取值服从泊松分布，其概率为：
[0065][0066]
采用伪随机数产生程序来确定每个10-kb的片段的k值。通过以上方法，为k562细胞的包含23条染色体的染色质区域构建了10个高分子链模型。图1-2以k562细胞的一个5mb区间(chr5:109-114mb)为例展示了这个区域的atac-seq的染色质可及性分布图和10个高分子链模型的珠子数目分布图。将这10个珠子数目分布进行平均得到的分布图与atac-seq的染色质可及性分布图是一致的。
[0067]
以上，为k562细胞的包含23条染色体的染色质区域构建了10个10-kb精度的携带染色质可及性信息的高分子链模型。注意，模型精度是指模型中染色质片段的长度。接下来，采用层级优化法模拟每个高分子链模型的随机折叠过程。过程如下：
[0068]
(1)在构建10-kb精度的高分子链模型时，先将染色质区域分割成连续的10-kb染色质片段，然后给每个10-kb片段分配珠子。现在，将相邻的两个10-kb染色质片段依次融合成20-kb的片段，融合后得到的20-kb片段的珠子数目等于融合前两个10-kb片段的珠子数目的平均数。这样，将10-kb精度的高分子链模型转变成20-kb精度的高分子链模型。依次类推，构建出了40-kb,80-kb,160-kb,320-kb,640-kb,1280-kb,2560-kb,and 5120-kb精度的高分子链模型。
[0069]
(2)基于5120-kb精度的高分子链模型产生一个随机结构，采用数值积分法模拟该随机结构的随机折叠。模拟中，随机结构的总势能包括三部分：采用的势能函数包含三部分：谐振子势能e
har
、排斥势能e
wca
和空间约束势能e
sphere
；其中，谐振子势能用于保证高分子链中相邻链珠之间的连续性，排斥势能用于避免两个不相邻的链珠之间发生重叠，空间约束势用于模拟由于细胞核的狭窄所产生的空间限制。这三类势能由多种差别不大的表达形式，本实施例采用以下形式：
[0070][0071][0072][0073]
其中r
n,n 1
是序列相邻的第n个和第(n 1)个两个珠之间的空间距离；r
n,m
序列不相邻的第n个和第m个两个珠子之间的空间距离；rn是第n个珠子的坐标；k是玻尔兹曼常数；t是温度；b是珠子的直径；σ＝b；n是聚合物模型中所有珠子的总数。模拟的步长为：
[0074][0075]
其中，f
max
是施加在所有珠子上的力中最大的力；m是珠子的质量。模拟中的每一步，按照玻尔兹曼分布随机产生每个珠子的速度。模拟200,000步后返回优化后的5120-kb精度3d结构。
[0076]
(3)在优化的5120-kb精度3d结构上，通过线性插值法，在相邻的珠子之间插入一个珠子，产生1个2560-kb精度的3d结构，以此结构作为初始结构进行第(2)步的优化，返回2560-kb精度的优化3d结构。以此类推，最后得到10-kb精度的优化3d结构。
[0077]
以上，构建了k562细胞的包含23条染色体的染色质区域的10个10-kb精度的携带染色质可及性信息的高分子链模型集合，并对集合中的每一个模型，都从50个不同的初始结构出发，进行50次层级优化，最后，得到包含500个10-kb精度的构象的集合，将其命名为hetero-ensemble。对hetero-ensemble中的每一个构象的每一对10-kb染色质片段质心之间的距离进行测量，构建出10-kb精度的单个距离矩阵(individual distance matrix)。设定片段之间的距离小于4个珠子直径的存在接触(contact)，将10-kb精度的单个距离矩阵转变为10-kb精度的单个接触矩阵(individual contact matrix)。将500个单个接触矩阵进行平均得到10-kb精度的平均接触矩阵(ensemble-average contact matrix)，该矩阵可以直接与多细胞的hi-c实验接触矩阵对照。用分离分数(separation score)分别识别出模拟的平均接触矩阵和hi-c接触矩阵中的tads，并进行比较。图3以k562细胞的一个5mb区间(chr5:109-114mb)为例，展示了hetero-ensemble的10-kb精度的模拟平均接触矩阵与实验hi-c接触矩阵的对比。可以看出这个区域中，hi-c实验检测到的绝大部分tad都被模拟矩阵预测了。图4还显示，在整个包含23条染色质的区域，实验检测到的tad边界有37％被模拟矩阵预测。图3和4表明，本发明能够在全基因组范围定量预测染色质tad结构。
[0078]
另外，用单细胞的fish成像数据来检验预测的hetero-ensemble。文献中报道了13,997个k562细胞的2-mb染色质区域(chr21:29.38-31.33mb)的30-kb精度的fish图像。基于这些实验数据，生成了该2-mb染色质区域的30-kb精度的fish平均距离矩阵。同时，根据模拟的10-kb精度的500个构象，计算了该区域的30-kb精度的平均距离矩阵。计算过程包括：(1)将10-kb精度构象中的相邻的3个10-kb染色质片段融合成30-kb的片段，并计算每个
30-kb片段质心的坐标；(2)测量每两个30-kb片段质心之间的距离，以此得到每个10-kb精度构象的30-kb精度的单个距离矩阵；(3)将500个单个距离矩阵进行平均，得到30-kb的平均距离矩阵。图5展示了30-kb精度的fish平均距离矩阵和30-kb模拟的平均距离矩阵。图6显示两者的相关系数高达0.935。再次证明，本发明能够预测染色质tad结构。
[0079]
进一步，从造血干细胞和祖细胞(hspc)到成熟免疫t细胞分化过程需要经历8个阶段，分别hspcs、mpp、clp、etp、、dn2、dn3、dn4和dp。以这八个阶段的细胞为例，根据dnase-seq实验数据，为这八类细胞的包含meis1基因的区域(chr11:15mb-25mb)构建了10条携带染色质可及性信息的高分子链模型。与生成hetero-ensemble集合完全一样的步骤，为每类细胞的该区域生成了一个包含500个10-kb精度的构象的集合。基于每个集合，计算了平均接触矩阵，并将它们与相应的多细胞hi-c接触矩阵对比。图7展示了计算的平均接触矩阵和hi-c实验接触矩阵。如图7所示，预测的接触矩阵显示，从dn2阶段过渡到dn3阶段时，染色体三维结构发生了显著的变化，是细胞命运的决定步骤。这一结论与hi-c接触矩阵显示的数据完全一致。这说明，本发明能够预测细胞分化过程中的关键染色质结构重组。