CN108220394B

CN108220394B - 基因调控性染色质相互作用的鉴定方法、系统及其应用

Info

Publication number: CN108220394B
Application number: CN201810011140.7A
Authority: CN
Inventors: 陈阳; 李炎剑; 贺毅; 张奇伟
Original assignee: Tsinghua University
Current assignee: Tsinghua University
Priority date: 2018-01-05
Filing date: 2018-01-05
Publication date: 2021-03-23
Anticipated expiration: 2038-01-05
Also published as: WO2019134586A1; CN108220394A

Abstract

本发明涉及一种运用多种基因组分析技术，例如RNA‑seq、ChIP‑seq或Hi‑C等手段，对经历状态转变(例如细胞分化)的样品的目标特性进行多方面表征和分析，从而识别出基因调控性染色质相互作用的方法，特别涉及一种鉴定在样品状态转变过程中，与样品状态转变相关的染色质相互作用和相关基因的方法。

Description

基因调控性染色质相互作用的鉴定方法、系统及其应用

技术领域

本发明涉及一种鉴定基因调控性染色质相互作用的方法，特别涉及一种鉴定在样品状态转变过程中，能够影响样品转变的染色质相互作用以及相应的效应基因的方法。

背景技术

染色质构象在基因表达的调控中起关键作用。研究发现有丝分裂间期染色体占据特定的领域，基因转录与基因相对于核纤层以及染色体领域的相对位置密切关联。最近使用高通量染色体构象捕获(Hi-C)的研究已经揭示，基因组被组织成几百个千碱基至一个兆碱基的拓扑相关结构域(Topologically Associating Domain，简称TAD)，并且TAD中的染色质区域更可能与同一个TAD内的其他区域发生作用，而不是与TAD以外的区域发生作用。并且在不同的细胞类型之间，大多数TAD位置是保持不变的，并显示出了进化上的保守性。

在同一个TAD的基因在面临激素刺激时或者在分化过程中显示出了相互协调的变化，这表明了TAD不仅是一个结构结构单元，而且还能够作为转录调节的功能单位。此外，在TAD内，由特定蛋白质或非编码RNA介导的长程染色质相互作用将远距离的调控区域，如增强子和基因启动子连接起来，从而使基因表达的远距离调控成为可能。

例如，在细胞分化中，常常伴随着关键基因的表达差异以及染色质三维结构或者构象的巨大改变，但是目前并没有有效的手段能够用于确定在这个过程中染色质结构的改变如何与关键基因的表达等行为的相互关联，以及这种相互关联是如何影响到细胞分化等状态变化的。因此本领域迫切需要一种新的方法，能够有效的分析并鉴定出在状态转变过程中，具有调节功能的染色质相互作用或者与受到染色质相互作用影响或调节的，对状态转变具有重要作用的关键基因或者调控因子。

发明内容

发明人经过长期的研究，获得了一种将染色质相互作用、特定基因的表达水平和/或识别位点与染色质内部基因调控这三个方面进行关联的方法，从而完成了本发明。

在第一个方面中，本发明涉及一种样品状态转变效应基因的鉴定方法，所述效应基因的表达受到样品状态转变中染色质相互作用改变的影响，其包括下列步骤：

(1)对处于第一状态和第二状态的样品进行比较，从而至少获得下列差异信息：基因可识别行为差异，以及存在于基因转录调控区域的染色质相互作用差异，和

(2)将步骤(1)获得的差异信息建立关联，获得与状态转变中转录调控区域的染色质相互作用差异有关的基因可识别行为差异，从而鉴定所述效应基因。

在一个实施方式中，其中所述样品是细胞。

在另一个实施方式中，其中所述基因可识别行为差异包括基因表达量的差异和/或基因调控区域基因组序列中结合模体分布的差异；优选的，所述基因表达量的差异是mRNA表达量差异或蛋白质表达量差异。

在另一个实施方式中，其中通过以下步骤获得步骤(1)中存在于基因转录调控区域的染色质相互作用差异：

(a)鉴定样品基因组中处于激活状态的启动子和/或增强子的位点；

(b)鉴定所有染色质相互作用的发生区域；

(c)整合步骤(a)和步骤(b)所获得的信息，得到位于激活状态的启动子和增强子之间存在的染色质相互作用，即存在于转录调控区域的染色质相互作用；和

(d)将不同样品之间存在于转录调控区域的染色质相互作用进行比较，得到存在于基因转录调控区域的染色质相互作用差异。

在另一个实施方式中，其中通过以下步骤获得步骤(1)中的基因可识别行为差异：

i)获得处于处于第一状态和第二状态的样品；

ii)取部分处于第一状态和第二状态的样品，分别进行转录表达分析，并比较样本间的mRNA表达量差异；优选的，转录表达分析采用RNA测序即RNA-seq法。

在另一个实施方式中，其中进一步包括：

iii)取部分处于第一状态和第二状态的样品，分别进行染色质开放区域序列的分析，优选使用ATAC-seq法，并分析染色质开放区域序列中所分布的转录因子结合模体，优选进一步比较处于第一状态和第二状态的样品间的结合模体的分布差异。

在另一个实施方式中，其中通过以下步骤分析步骤(1)中存在于基因转录调控区域的染色质相互作用差异：

iv)取部分处于第一状态和第二状态的样品，鉴定分别处于激活状态的启动子和/或增强子信息，优选的，采用ChIP-seq法进行鉴定，所述ChIP-seq法中所使用的抗体优选为H3K4me3和H3K27ac的结合抗体，所述抗体分别结合H3K4me3和H3K27ac形成信号峰，分别代表了处于激活状态的启动子和增强子位点；

v)另取部分处于第一状态和第二状态的样品，采用染色质构象捕获技术，优选采用高通量染色质构象捕获技术，例如Hi-C法、insituHi-C法、BL-Hi-C法或ChIA-PET法，获得全基因组染色质相互作用的信息；或者利用4C或5C法获得局部染色质相互作用的信息；

vi)将参考基因组序列划分成一定大小的区域，优选的，所述区域大小在1-40kb之间，例如1kb、5kb、10kb、15kb、20kb、25kb、30kb、35kb或40kb，基于步骤iv)获得的活性状态的启动子和增强子信息，通过比对分别获得包含有活性的启动子和增强子位点的区域，优选将包含启动子的区域命名为基因区域，将包含增强子序列的区域命名为调控区域；

随后，结合步骤v)获得的染色质相互作用的频率信号，识别出发生于基因区域和调节区域之间的染色质相互作用频率信号，从而得到基因调控相关的染色质相互作用；然后将处于第一状态和第二状态的样品之间的基因调控相关染色质相互作用进行比较，其中具有统计显著性差异的，被鉴定为基因调控性染色质相互作用差异(differential gene-regulatory interaction)，包括相对于处于第一状态的样品，处于第二状态的样品中增强的基因调控性染色体相互作用和/或减弱的基因调控性染色体相互作用。

在另一个实施方式中，其中步骤(2)具体包括下述步骤：

a)将基因表达量差异与基因调控性染色质相互作用差异结合，选择在不同状态的样品中，存在于增强或减弱的基因调控性染色质相互作用内部、同时表达量也具有显著变化的基因；或

b)将基因组转录调控区域中转录因子结合模体分布的差异与基因调控性染色质相互作用差异信息相结合，选择在不同状态的样品中，存在于增强的或减弱的基因调控性染色质相互作用内部、转录因子结合模体分布也发生显著变化的基因；或

c)将基因表达量差异、转录因子结合模体分布的差异与基因调控性染色质相互作用差异信息相结合，选择在不同状态的样品中，存在于增强的或减弱的基因调控性染色质相互作用内部、基因组转录调控区域中结合模体分布发生显著变化、同时表达量也具有显著变化的基因。

在另一个实施方式中，此外还包括对筛选出来的效应基因进行功能研究以确定其功能的步骤。

在另一个实施方式中，其还进一步包括鉴定获得基因调控性染色质相互作用的步骤，即将能够影响步骤(2)鉴定得到的效应基因表达的染色质相互作用，作为基因调控性染色质相互作用。

在另一个实施方式中，其中所述样品状态转变通过下述方式实现：化学试剂诱导、自然分化和/或物理刺激。

在第二个方面中，本发明涉及一种鉴定能够调控样品状态转变的染色质相互作用的方法，其包括第一方面中任一实施方式所述的步骤。

在第三个方面中，本发明涉及一种鉴定参与样品状态转变所涉及染色质相互作用的调控因子的方法，其包括第一方面中任一实施方式所述的步骤。

在第四个方面中，本发明涉及一种鉴定能够调控染色质相互作用的物质的方法，其包括：利用第一个方面的任一实施方式鉴定得到染色质相互作用的效应基因或基因调控性染色质相互作用，随后将待测物质与样品接触，分析所述效应基因或基因调控性相互作用的变化。

在第五个方面中，本发明涉及一种样品状态转变效应基因的鉴定系统，所述效应基因的表达受到样品状态转变中染色质相互作用改变的影响，包括下述模块：

(1)基因可识别行为差异分析模块；

(2)转录调控区域的染色质相互作用差异的分析模块；和

(3)效应基因鉴定模块；

所述系统能够获得与染色质相互作用差异相关的基因可识别行为差异，从而获得受染色质相互作用影响的样品状态转变效应基因。

优选的，所述系统还进一步包括基因调控性染色质相互作用鉴定模块，从而鉴定能够影响所述效应基因表达的基因调控性染色质相互作用。

在一个实施方式中，其中所述基因可识别行为差异分析模块能够分析基因表达量的差异和/或基因的转录调控区域的基因组序列中转录因子结合模体分布的差异。

在另一个实施方式中，其中转录调控区域的染色质相互作用差异的分析模块能够执行下述分析：

(b)鉴定所有染色质相互作用的发生区域；

(d)将不同样品之间存在于转录调控区域的染色质相互作用进行比较，得到存在于转录调控区域的染色质相互作用差异。

在又一个实施方式中，所述效应基因鉴定模块能够执行下述分析：

在第六个方面中，本发明涉及一种检测试剂盒，其包含本发明第一至第四方面所述方法中所使用的试剂。

本发明的方法通过整合多个组学试验方法和结果，从而建立了染色质构象与转录调控之间的分析和鉴定方法。该方法可应用于分析多个生物学过程，如细胞分化、个体发育、细胞变异、疾病治疗等方面，从而可以在染色质构象水平鉴定出对上述过程具有重要影响的染色质相互作用和调控因子。

附图说明

图1显示了本发明一个实施方式的整体流程。

图2显示了本发明一个实施例中，经过全反式维甲酸(All-trans-retinoic acid，缩写为ATRA)诱导后的HL-60细胞与对照组相比mRNA表达变化，图2A显示了维甲酸诱导后基因表达的变化情况，图2B显示了差异表达基因在不同GO分类中的富集情况。

图3显示了经过ATRA诱导后，TAD内部的染色质相互作用频率变化(图3A)以及频率变化与差异表达基因的关系(图3B)。

图4显示了具有TAD内染色质相互作用频率变化与H3K4me3和H3K27ac修饰变化的关系。

图5显示了本发明一个实施例中所获得的基因表达调控相关的染色质相互作用差异，其中图5A显示了H3K4me3信号和H3K27ac信号在ATRA处理组和对照组中的变化，图5B显示了ATRA和对照组的特异性峰和共有峰分别与转录起始位点的距离，图5C显示了基因调控相关的染色质相互作用的寻找方式，图5D显示了Gain与Loss组的基因区域内H3K27ac的信号，以及在基因区域内H3K4me3的信号的相对比较，图5E显示了Gain和Loss组中差异表达基因的情况。

图6显示了本发明一个实施例中经过ATRA诱导后，ATAC-seq测定染色质开放区域的结果比较。图6A和6B显示了各组的特异性峰以及峰在基因组的分布，图6C和图6D显示了分别在调控区域和基因区域中，Gain组和Loss组相互作用内ATRA处理组和对照组特异性以及共有的染色质开放信号分布，其中Gain组更富集ATRA处理组特异性信号，而Loss组更富集对照组特异性信号，从中还能够看出调控区域相比于基因区域更富集了特异性信号。图6E显示了在对照组和ATRA诱导组中，在染色质开放区域上富集的转录因子结合模体(motif)。

图7显示了在Gain组(图7A)和Loss组(图7B)中的染色质开放区域内具有结合模体的转录因子(TFs)分布比较，从中可以发现，GATA结合模体在两组间具有特异性差别。

图8A显示了通过转录因子-靶基因调控网络分析，表明GATA2位于网络核心枢纽，图8B显示了ATRA诱导HL-60细胞分化进程中，GATA2与其他富含ATAC-seq峰的转录因子之间的相互作用关系。

图9A和图9B分别显示了ATRA诱导后GATA2基因区域和调节区域相互作用和GATA2基因表达变化；图9C和图9D分别显示了ARTA诱导后ZBTB16基因区域和调节区域相互作用和ZBTB16基因表达变化；

图10A显示了GATA2基因上游4C染色质构象捕获的实验结果，图10B显示了验证染色质环结构的原位FISH结果，图10C显示了FISH实验中红绿两种荧光分布的皮尔森相关系数。

图11A显示了ZBTB16基因区域4C染色质构象捕获的实验结果，图11B和图11C分别显示了Gata2和Zbtb16在诱导分化中的模型，其中解释了染色质结构、转录因子结合和基因表达的关系。

具体实施方式

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

定义

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

术语“状态转变”是指对于同一样品而言，通过某种特定的外加诱导或者内部的自然过程而导致样品的性质或形态等发生改变。例如化学试剂的诱导分化、物理刺激或细胞在自然生理过程中的分化，如响应外部激素或其他信号分子或者细胞内部基因或蛋白的作用，而引发的细胞自然分化过程。在一个实施方式中，本发明中“至少两种处于不同状态的样品”就是通过状态转变而形成的。

“处于第一状态的样品”和“处于第二状态的样品”是指经过状态转变过程后得到的两种不同状态的样品。其中在一些实施方式中，“处于第一状态的样品”是状态转变前的样品，而“处于第二状态的样品”是经过状态转变后的样品。

术语“效应基因”是指参与到样品状态转变过程的基因，效应基因可能是样品发生状态转变的原因，即该基因能够引发状态转变的过程，例如在细胞分化诱导模型中，该基因可以是直接响应外加诱导的基因，从而引发细胞的分化；另外该基因也可以是状态转变过程中间的环节，或者仅仅作为状态转变的结果。需要注意的是，在本文中“基因”可以指基因、也可以指基因的表达产物，如mRNA转录本或者蛋白质，在一个实施方式中，效应基因可以是转录因子。

术语“转录调控区域”，或者叫做调控区域，是指对基因组DNA中位于基因的上游或下游一定的范围内，例如10kb-1Mb，50kb-500kb，100kb-200kb范围内，包含启动子、增强子等反式作用因子(例如转录因子)结合位点的一段区域。

在本发明一个实施方式中，转录调控区域的分析对从大量的可供选择的候选基因中挑选出目标效应基因具有重要作用，例如确定基因的转录调控区域内是否存在相互作用的启动子和增强子，以及转录调控区域是否是开放序列，以及存在哪些转录因子的结合模体；在不同状态的样品中，上述增强子和启动子之间的相互作用是否发生改变，同时结合模体的分布是否也发生了改变等，再与其他信息结合(如基因表达量)就可以有效的获得效应基因。

术语“基因调控相关的染色质相互作用”，也叫做“存在于转录调控区域的染色质相互作用”，是指在基因转录调控区域的不同调控元件内，例如启动子和增强子的序列之间，发生的染色质相互作用。

在本发明的一个实施方式中，将参考基因组序列划分成一定大小的区域，所述区域的大小可以根据染色质构象分析的数据深度例如测序深度进行调整，优选的范围是1kb-40kb，例如1kb、2kb、3kb、4kb、5kb、6kb、7kb、8kb、9kb、10kb、11kb、12kb、13kb、14kb、15kb、16kb、17kb、18kb、19kb、20kb、21kb、22kb、23kb、24kb、25kb、26kb、27kb、28kb、29kb、30kb、31kb、32kb、33kb、34kb、35kb、36kb、37kb、38kb、39kb或40kb。在一个具体实施方式中，所述区域大小为40kb。

接下来，基于前述步骤iv)获得的活性状态的启动子和增强子信息，通过比对分别获得包含有活性的启动子和增强子位点的区域，优选将包含启动子的区域命名为基因区域，将包含增强子序列的区域命名为调控区域；随后，结合步骤v)中的染色质相互作用的信号，分析基因区域和调节区域之间的染色质相互作用频率信号(即特定区域之间的染色质相互作用数，例如具体可表现为Hi-C数据中两端分别落在特定区域的读段数)，当基因区域和调节区域之间具有可识别强度的接触信号时，被认为存在基因调控相关的染色质相互作用。

术语“基因调控性染色质相互作用差异”，在本文中有时候也被称为“转录调控区域的染色质相互作用差异”、“位于转录调控区域的染色质相互作用差异”或“存在于转录调控区域的染色质相互作用差异”，其表达的含义相同。其获得方法是：比较处于第一状态和第二状态的样品的基因调控相关的染色质相互作用，其中具有显著性差异的，被鉴定为基因调控性染色质相互作用差异(differentialgene-regulatory interaction)，其中的“显著性”优选是指统计学意义上的显著。例如当采用假设检验时，p<0.05或p<0.01。

术语“基因调控性染色质相互作用”，是指在处于第一状态的样品和处于第二状态的样品间具有显著性差异的那些基因调控相关的染色质相互作用，即属于染色调控相关的染色质相互作用的一部分，在鉴定出的基因调控性染色质相互作用中，染色质相互作用与基因调控之间的关联更加确定。

实际上，前述鉴定得到的“基因调控性染色质相互作用差异”就可以认为是“基因调控性染色质相互作用”的比较结果。根据所述“差异”的不同类型，基因调控性染色质相互作用也分为两种类型：即处于第二状态的样品相对于处于第一状态的样品增强的基因调控性相互作用(在本发明的实施例中，这一类相互作用也被归为Gain组)和/或减弱的基因调控性相互作用(在本发明的实施例中，这一类相互作用也被归为Loss组)。

术语“基因可识别行为差异”指在不同状态的样品中与基因的性质、状态等相关的可以被定性或定量观测到的差异。而其中所述“基因”并非特定的某一部分或经过人工干预而预先选定的基因，而是在定量或定性的分析中，被观察到具有可识别行为差异的整体基因集合。有时候，为了清楚，可以将上述“基因”定义为备选基因或候选基因，但需要注意的是，在本文中即使使用“备选基因”或“候选基因”的表述，也并不表明这些备选基因或候选基因是需要人为预先选定范围的一部分。

术语“结合模体”是指存在于基因组DNA上，能够被转录因子等反式作用因子结合从而调节目标基因，例如调节本发明中效应基因表达的元件。

术语“结合模体分布的差异”是指：在不同状态的样品之间，整体或部分结合模体在数量、位置或存在与否方面的差异，或者，指位于感兴趣的区域中的部分或特定结合模体在数量、位置以及存在与否方面的差异。

术语“染色质相互作用”是指染色质不同的位点之间的长距离相互作用，从而形成染色质的高级构象从而维持染色质结构或者促进基因的表达。

术语“染色质相互作用频率”，在本文中也叫做“Hi-C相互作用频率”或Hi-Ccontact frequency，是指在进行染色质构象分析时，寻找染色质相互作用时发现的不同区域之间相互作用的信号，表现为Hi-C数据中两端分别落在特定区域的读段数。

术语“染色质开放区域序列”是指染色质中由于无核小体结合等原因而暴露的DNA序列，并可被反式作用因子如转录因子结合。

术语“ChIP-seq”是指将免疫共沉淀技术(ChIP)与高通量测序相结合的技术，从而高效地在全基因组范围内检测与组蛋白、转录因子等相互作用的DNA区段。其原理在于：首先通过染色质免疫共沉淀技术(ChIP)特异性地富集目的蛋白结合的DNA片段，并对其进行纯化与文库构建；然后对富集得到的DNA片段进行高通量测序。随后将获得的数百万条序列读段精确定位到基因组上，即获得全基因组范围内与组蛋白、转录因子等互作的DNA区段信息。

术语“染色质构象捕获技术”是指所有能够实现染色质不同空间位置之间的关系从而建立染色质三维结构信息的技术，其包括普通的3C技术，即Chromosome ConformationCapture，也包括结合高通量测序的染色质构象捕获技术。

术语“高通量染色质构象捕获技术”，是指结合高通量测序技术以及生物信息分析方法，从而有效分析全基因组范围内整个染色质DNA在空间位置上的关系，并获得高分辨率的染色质三维结构和染色质相互作用信息的方法。在本文中，该技术至少包括Hi-C、基于Hi-C的改进型技术in situ Hi-C，以及在in situ Hi-C法的基础上进一步引进bridge-linker后获得的BL-Hi-C法，另外ChIA-PET法在本文中也属于一种高通量染色质构象捕获技术。

术语“ATAC-seq”是指一种分子生物学中研究染色质可接近性的技术，由ATAC实验和高通量测序两部分组成，ATAC-seq实验的关键部分是转座酶Tn5对样品基因组DNA的作用。转座子优先并入一般没有核小体(无核小体区域)或暴露DNA段的基因组区域。因此，基因组中某些基因座序列的富集表明该区域不存在核小体，处于DNA结合蛋白等核机器能进入的松散暴露状态，提供有关染色质区段转录活跃状态的信息。ATAC-seq采用突变的多活性转座酶，允许高效切割暴露的DNA和同时连接特定序列的接头。分离接头连接的DNA片段，通过PCR扩增后用于高通量测序。

实施例

下述实施例以ATRA诱导HL-60分化为例，示例性的说明本公开方法如何寻找在上述分化过程中具有重要调节的基因并进行相应的分析。需要说明的是，本领域技术人员可以理解，本公开中的方法并不限于所述实施例所示例的方法，而是适用于任何处于两种不同状态下的样本中，相关目标调节基因的寻找和分析。

实施例1细胞培养和ATRA诱导

HL60细胞购自国家实验细胞资源共享平台(中国北京协和医学院)。细胞维持在补充有10％胎牛血清(FBS，Gibco，USA)，50单位/mL青霉素和链霉素(Gibco，USA)和非必需氨基酸(Gibco，USA)的RPMI-1640培养基(Gibco，美国)中。

对于粒细胞分化，2x10⁵/ml HL-60细胞使用1μM ATRA(溶于乙醇的1mM储备液，Sigma，USA)处理4天(称为ATRA组)；用等量乙醇处理的细胞被称为对照组。第2天更换培养基，同时加入ATRA/乙醇。

实施例2获得对照组和ATRA组细胞的RNA差异表达信息

方法：RNA seq

过程

利用TRIZOL(Ambion，USA)法从对照细胞以及经过ATRA处理的细胞中提取总RNA。文库构建和测序均由安诺优达公司(中国)进行。

对于RNA-seq分析，首先去除接头序列(adapter)，随后用Bowtie将数据比对回参考基因组hg19，滤掉核糖体RNA的测序读段(reads)。经过上述步骤后，将剩余的读段数据与转录组数据RSEM v1.2.7进行比对，并进行定量分析。注释文件下载于加利福尼亚大学圣克鲁斯分校(UCSC)基因组浏览器(genome browser)的人类基因组hg19版本(humanhg19assembly)。使用Deseq2 1.4.5版本的软件包，根据基因方差分布估计(gene-wisedispersion estimates)的平均值计算差异基因表达。基于调整p值等于0.01和倍数变化值取log2后大于0.9，确定表达具有显著差异的基因。Gene Ontology分析使用DAVID。使用R3.3.1版本的ggplot2软件包绘制RNA-seq图(plots)。

结果

差异表达基因(DEGs)分析显示在ATRA诱导后有941个上调基因和611个下调基因(图2A)。GO分析显示，“免疫反应”和“白细胞活化”类基因被显著富集，这实际上与嗜中性粒细胞的终末分化过程一致(图2B)。

实施例3TAD内部染色质相互作用水平与基因表达水平、基因启动子和增强子活性密切相关

方法：BL-Hi-C

过程

文库构建：用1％甲醛处理细胞以交联细胞中的蛋白质以及蛋白质和DNA，然后使用裂解缓冲液(50mM HEPES-KOH，150mM NaCl，1mM EDTA，1％Triton X-100和0.1％SDS)重悬。然后用酶HaeIII将基因组消化成具有平末端的片段。将DNA片段平末端用腺嘌呤处理，并在16℃下与含有生物素的桥连接物(bridge linker)连接4小时，并用外切核酸酶(NEB)消化未连接的DNA片段。接下来，细胞用蛋白酶K(Ambion)消化过夜，利用酚-氯仿(Solarbio)结合乙醇沉淀提取并纯化DNA。然后，使用S220聚焦超声波仪(Covaris)将DNA片段化，并通过链霉亲和素包被的Dynabeads M280(Thermo Fisher)结合生物素标记的DNA片段。对由磁珠制备的文库进行Illumina测序并通过PCR进行扩增。用AMPure XP珠(Beckman，德国)纯化后，使用Illumina HiSeq 2500测序仪进行测序。

Hi-C数据分析：首先，将读段中的桥连接物序列(序列：CGCGATATCTTATCTGACT或GTCAGATAAGATATCGCGT)去除，如果完整接头将读段分成两个片段，则保留5'片段。其次，经过上述处理的读段被对应到人类基因组hg19版本，同时去除重复片段。第三，基于DNA链信息估计读段对的距离阈值，结合读段对和链的信息，将相互作用对分成下面几种类型：完整片段(读段对内部无发生连接发生)、自连接、染色体间连接和染色体内连接。对于识别4碱基对的限制性内切酶而言，染色体内连接的阈值是大约3kb。

Hi-C数据校正：使用迭代校正法(ICE)来纠正系统偏差，之后产生40kb分辨率的相互作用矩阵。

对于TAD分类，我们使用基于聚类的Hi-C结构域搜索(CHDF)法进行。

计算TAD相互作用倍数改变：为了计算ATRA诱导时TAD内部和外部Hi-C计数的倍数变化，首先计算重复之间的倍数变化。对于每个TAD，将对照组和ATRA组细胞的倍数变化相结合以产生背景分布(分别计算内部和外部倍数变化)。然后，将ATRA处理细胞与对照细胞之间的倍数变化引入背景分布中以获得基于它们在背景分布中位置的p值。在两个重复中p值<0.05的TAD被定义为显着变化的TAD。

方法：ChIP-seq文库构建和数据分析

过程

利用1％甲醛对经过ATRA处理的HL-60细胞和对照细胞进行交联处理。然后，使用裂解缓冲液(50mM HEPES-KOH，150mM NaCl，1mM EDTA，1％Triton X-100,0.1％脱氧胆酸钠和1％SDS)裂解细胞膜。染色质用FA裂解缓冲液(50mM HEPES-KOH，150mM NaCl，1mM EDTA，1％Triton X-100，0.1％脱氧胆酸钠和0.1％SDS)重悬，随后利用超声处理处理器(Cole-Parmer，美国)进行片段化。利用与H3K4me3和H3K27ac抗体(Abcam，England)预温育的Dynabeads(Thermo Fisher)过夜处理进行免疫沉淀。洗涤和纯化DNA后，使用TruePrep DNA文库制备试剂盒(Vazyme，中国)根据制造商提供的实验手册进行文库构建。使用IlluminaHiSeq 2500测序仪对文库进行测序。

对于ChIP-Seq分析，首先去除接头序列。然后，使用Bowtie将测序读段比对到人类基因组hg19版本上。组蛋白修饰峰使用MACS2软件调用生成，参数设置为'-g hs--nomodel--broad'。在两次重复中均存在的峰(bedtools软件，1bp最小重叠)被认为是置信峰。然后，比较对照和ATRA处理的细胞所获得的置信峰，从而区分ATRA处理组特异峰、对照组特异峰和重叠峰。峰值比较使用bedtools中的intersectBed软件进行。

分析

根据TAD内部Hi-C相互作用频率变化对包含表达基因(n＝3362)的TAD进行排序，并根据Hi-C相互作用频率的统计分布将这些TAD分类，并计算这些TAD中表达基因的分布情况(图3A，3B)。结果表明，经过ATRA诱导后，相比于内部相互作用变化并不明显的TAD，内部相互作用发生明显变化的TAD中更可能具有差异表达的基因。而且，内部染色质相互作用增加的TAD会富集上调的差异表达基因，而内部染色质相互作用减弱的TAD则表现出富集下调的差异表达基因(图3B)。以上结果显示，TAD内基因表达变化与内部染色质相互作用频率之间存在正相关性。

进一步的，为了定性表观表征TAD内的表观遗传状态，我们使用针对H3K4me3和H3K27ac的抗体分别在对照和ATRA处理的细胞中进行ChIP-seq，其主要标记活性启动子和增强子。通过计算TAD内的ChIP-seq信号变化，可以发现内部染色质相互作用增加的TAD的H3K4me3和H3K27ac的水平也增加，而染色质相互作用减少的TAD则显示相反的变化(图4)。这提示了TAD内部表观遗传学活性的改变可能会影响染色质构象变异和差异基因表达。

通过本实施例可知，TAD内部染色质相互作用水平与基因调控非常密切，可以进一步用于寻找目标功能基因。

实施例4“基因表达调节相关的染色质相互作用差异(differential gene-regulatory chromatin interactions)”可有效指示基因表达差异

在实施例3的基础上，为了更好的将TAD内部染色质相互作用水平定量化从而应用于基因的鉴定，首先利用实施例3中的ChIP-seq数据，并在ATRA组细胞中鉴定了12295个H3K4me3峰和12493个H3K27ac峰，在对照组细胞中鉴定了14263个H3K4me3和22149个H3K27ac峰(图5A)。其中H3K4me3峰代表活性的启动子，H3K27ac峰代表活性的转录区域与增强子。通过计算每个H3K27ac峰的最接近转录起始位点(TSS)的距离，我们发现ATRA组和/或对照组特异性峰，相较于二者的共有峰，常常更加远离TSS(图5B)，表明在ATRA诱导后出现更远端的调节变化。

染色质相互作用能够将远端的调节元件、启动子和转录起始位点拉近从而用于转录起始，为了更进一步确定究竟染色质相互作用的改变如何影响了基因表达以及影响了哪些基因表达，进行了如下的步骤：

首先将整个基因组分成大小为40kb的bin。接下来，如果bin包含表达基因的启动子(即H3K4me3峰)，该bin被标注为“基因区域”；如果包含启动子远端的H3K27ac峰(图5C)，则该bin被标注为“调控区域”。然后，利用基因区域和调节区域之间的Hi-C相互作用关系(Hi-C读段)来表示基因调控相关的染色质相互作用的强度。由此，我们使用Hi-C数据以40kb的分辨率生成了染色质相互作用矩阵，并且矩阵的相应元素中的计数表示基因和调控区域之间的相互作用强度。例如，如果第i个bin是基因区域，第j个bin是调控区域，则[i，j]的计数代表基因调控相关的染色质相互作用。相同TADs内的基因调控相关的染色质相互作用作为差异相互作用分析的输入。差异相互作用分析是基于MA曲线法和随机抽样模型而进行的，即：

认为Hi-C实验是对多数细胞中染色质相互作用的采样，因此Hi-C数据中两区域相互作用的读段遵循二项分布。使C₁和C₂分别表示从对照和ATRA处理的细胞中得到的特定基因调控性染色质相互作用的计数，其具有C_i～二项式分布(n_i，p_i)，i＝1,2，其中n_i表示Hi-C计数的总数，p_i表示来自该基因调控性染色质相互作用的计数的概率。我们定义M＝(log₂C₁-log₂C₂)/2，并且A＝(log₂C₁+log₂C₂)/2。在随机抽样的假设下，假设A＝a(a是A的一个观察值)，M的条件分布遵循近似的正态分布。对于MA图上的每个基因调控性染色质相互作用，我们进行H₀：p₁＝p₂与H₁：p₁≠p₂的假设检验。然后，基于条件正态分布分配p值。使用R3.3.1版本的DEGseq软件包进行分析，参数设置为“MARS”。

结果：根据在对照和ATRA处理的细胞之间Hi-C接触差异是否显著(基于Benjamini校正的p<0.001)，挑选出了基因调控性染色质相互作用差异，并鉴定了422对增强的基因调节相关的染色质相互作用(可简称为“Gain组”)和330对减弱(或降低)的基因调节相关的染色质相互作用(可简称为“Loss”组)(图5C)。

进一步的，在Gain组中，相应基因区域中H3K27ac信号的倍数变化与相应基因区域中的H3K4me3信号的倍数变化均显著增加，而在Loss组中上述两组信号则呈现相反趋势(图5D)。这些结果可以至少部分地解释活性组蛋白状态和染色质相互作用强度之间的正相关性以及与差异基因表达的正相关性。经统计，分别共计有430个和323个基因分别涉及了增强的基因调节相关的染色质相互作用以及降低的基因调节相关的染色质相互作用；进一步与RNA表达数据结合可以发现，当ATRA诱导时，上述基因中分别有164个和61个基因显示出了差异性的表达(图5E)，其比例显著高于整体比例(在12266个基因中具有1649个差异表达基因)可见，Gain组或Loss组的相互作用分别在很大程度上能够导致涉及的相应基因的上调或下调。

实施例5基因调节性染色质相互作用差异与染色质的可接近性之间存在关联

方法：ATAC-seq

流程

将对照样品和ATRA处理的样品分别在裂解缓冲液(10mM Tris-HCl(pH 7.4)，10mMNaCl，3mM MgCl₂和NP-40)中置于冰上10分钟以制备细胞核，在细胞裂解后立即旋转以除去上清液。随后将细胞核与Tn5转座体和标记缓冲液在37℃下孵育30分钟(Vazyme，中国)。在标记后，直接将停止缓冲液加入到反应体系中以结束标记反应。随后进行12个循环的PCR以扩增文库。PCR反应后，使用1.2×AMP珠(Beckman，德国)纯化文库，然后使用IlluminaHiSeq 2500测序仪对得到的文库进行测序。

去除接头序列并且将序列比对到人类基因组hg19版本中。利用MACS2使用默认参数调出ATAC-seq峰值，随后进行峰值比较。为了鉴定富含ATAC-seq峰的序列基序，使用HOMER程序中的MififsGenome.pl。AnnotatePeaks.pl用于识别包含某些基序的特定峰。ATAC-seq峰的GREAT分析参见文献(McLean CY,Bristor D,Hiller M,Clarke SL,SchaarBT,Lowe CB,et al.GREAT improves functional interpretation of cis-regulatoryregions.Nat Biotechnol.2010May2；28(5):495–501)。

结果

计算位于具有差异性基因调节性染色质相互作用的基因区域中和调节区域中的ATAC-seq峰，发现增强(Gain)和减弱(Loss)的相互作用中分别富集了ATRA特异性和对照特异性峰(图6C和6D)。调控区域显示出比基因区域更强的富集倾向，表明在开放染色质区域，特别是远端调控区域中TF结合的变化，调节了染色质相互作用的形成。

为了表征ATRA诱导后的转录因子结合状态，我们使用HOMER软件针对对照组和ATRA诱导组的特异性ATAC-seq峰进行基序分析(具体操作见实施例5)。包括CTCF，PU.1，RUNX和CEBP的大多数转录因子在ATRA处理的细胞或对照细胞之间高度相似(图6E)。通过ATAC数据分析观察到PU.1mRNA表达具有轻微上调(～1.9倍)；在RUNX家族成员中，只有RUNX3在ATRA处理后显示出显著的mRNA水平变化(～3.2倍，上调)，表明其可能在ATRA诱导中起调节作用。值得注意的是，GATA结合基序仅在对照细胞中富集(图6E)，并且GATA2mRNA的表达水平在分化后显著下调(～0.06倍)，这表明失去GATA2的结合与ATRA诱导过程相关。

实施例6在ATRA诱导模型中，通过整合多组学数据鉴定与染色质结构改变相关的关键基因转录

Hi-C和ChIP-seq的实验操作和结果分析具体见实施例3和4，从而获得增强的基因调节相关的染色质相互作用组(Gain)(图7A)或减弱的的基因调控性染色质相互作用组(Loss)(图7B)。然后在上述两个相互作用组中，对ATAC-seq峰利用HOMER软件进行motif分析(具体见实施例5)。

结果显示，具有结合位点的PU.1、RUNX JUNB等转录因子均在两组中存在富集现象，但在Gain组和Loss组间不存在富集有显著差异的现象，这表明这些转录因子与改变染色质相互作用的关联较低。但值得注意的是，GATA模体序列(如GATA1和GATA2)只在Loss组的调节区域中显著富集，这表明了GATA转录因子在染色质相互作用中具有独特的作用，进一步的结合RNA-seq的数据，由于GATA2 mRNA的表达水平在分化后显著下调(～0.06倍)，GATA2结合的丧失可能有助于ATRA诱导过程；由此，通过综合上述多组学的数据，成功鉴定得到GATA2作为候选基因。

实施例7更多相关转录因子的发现

为了进一步描述转录因子和差异表达基因在分化过程中的关系，我们将上述转录因子和差异表达基因定位到HTRI TF-Target网络(图8A)。在子网络中，我们发现GATA2是网络程度最高的枢纽节点。此外，GATA2与大多数富含ATAC-seq峰的转录因子和已知的粒细胞分化的关键调控因子(图8B)显示了相互作用。总之，通过整合染色质可及性信息和转录调控网络，我们发现GATA2可能在ATRA诱导的HL-60分化过程中起到重要的转录因子的作用。

实施例8ATRA诱导降低了GATA2启动子和调控区域之间的染色质相互作用

根据差异基因调控染色质相互作用分析，我们观察到在ATRA刺激后含有Gata2的基因区域和上游调节区域之间的相互作用显著降低(图9A和图9B)。为了详细描述与Gata2相关的染色质构象，接下来使用了Gata2启动子作为诱饵，在对照和ATRA处理的细胞中进行了4C染色质构象捕获技术。

具体流程如下：

首先，用1％甲醛交联细胞，并通过重悬于裂解缓冲液(500μl10mMTris-HClpH8.0,10mM NaCl，0.2％Igepal CA-630和50μl蛋白酶抑制剂)中分离细胞核，然后，用1×NEBuffer2(NEB，英国)洗涤细胞核，并于65℃下使用0.3％SDS处理。随后利用HindIII于37℃酶切过夜，然后在4℃进行近端连接4小时。用蛋白酶K(Ambion，USA)处理反交联后，使用酚-氯仿(Solarbio)提取结合乙醇沉淀从而纯化DNA。第二次酶切步骤使用DpnII过夜进行。然后，我们进行第二次连接和乙醇沉淀以提取DNA。最终使用QIAquick PCR纯化试剂盒(QIAGEN，德国)根据制造商的方案纯化DNA。PCR反应后，我们使用AMPure珠(Beckman，Germany)纯化4C文库，并使用Illumina HiSeq 2500测序仪对文库进行测序。

对于4C-seq数据分析，首先使用SAMtools中的cutadapt软件移除衔接子序列。使用Bowtie将读长对应到人类hg19版基因组信息。然后，使用RPM标准化，使用R3.3.1中的r3Cseq包处理对应后的数据。

结果

在对照细胞中，我们发现Gata2启动子与3个被H3K27ac峰所指示的上游增强子(即chr3：128240590-128254410,128262419-128292429和128309790-128334446)具有强的染色质相互作用。增强子E3(GATA2启动子上游大约80kb)的位置非常接近已知的增强子，证实了4C数据的可靠性。在ATRA诱导后，Gata2启动子和上游区域之间的相互作用强度在不同程度上下降(分别是E1，E2和E3的0.54-，0.46-和0.4倍)，与H3K27ac的下降一致(图10A)。模体分析显示这些区域分别含有PU.1/RUNX1和GATA模体序列，并且仅在对照细胞中观察到开放染色质状态(图10A)。关键转录因子(PU.1/RUNX1/GATA)在远端调控区域结合的缺失破坏了染色质环并抑制了GATA2表达。为了进一步证实ATRA诱导后染色质环的消失，我们在对照和ATRA处理的细胞中进行了三维DNA荧光原位杂交(FISH)，结果表明对照细胞相对于ATRA处理的细胞，结合Gata2的启动子和增强子E3的探针信号具有更大的重叠(图10B)。此外，两种探针之间的Pearson相关系数在对照细胞中高于ATRA处理的细胞中，进一步验证ATRA诱导后导致的染色质环的破坏(图10C)。

实施例9染色质相互作用的变化(失去染色质环)抑制了ZBTB16 mRNA的表达

在此前提到的转录调控网络中(图8A)，另一个编码锌指蛋白ZBTB16(也称为PLZF)的重要基因涉及了基因调节性染色质相互作用差异。在ATRA诱导后，该基因的表达量和染色质相互作用显着降低(图9C)。之前的一项研究(Tang Z,Luo OJ,Li X,Zheng M,Zhu JJ,Szalaj P,et al.CTCF-Mediated Human 3D Genome Architecture Reveals ChromatinTopology for Transcription.Cell.2015Dec17；163(7):1611–1627)使用ChIP-PET在K562细胞中鉴定了Zbtb16基因位点附近被CTCF结合的3个上游和2个下游锚点(anchor)。为了确定在ATRA诱导下ChIA-PET发现的染色质环是否改变，我们使用5种ChIA-PET锚作为诱饵进行4C分析(图11A)。其中一个3'锚点的结果显示Zbtb16的5'和3'之间的染色质环消失，同时伴随着ZBTB16mRNA的显着降低(图9D)。这个结果与先前预测相符，即染色质的5'和3'环保持高基因表达水平。在ATRA处理和对照细胞中，在Zbtb16的5'处观察到ATAC-seq峰。出乎意料的是，只有ATRA处理的细胞中的峰富集了PU.1基序，表明PU.1可以在ATRA处理后与5'锚结合。

基于实施例8和9的结果，我们提出了两个模型来解释在GATA2和Zbtb16区域中ATRA诱导分化过程中的染色质结构、转录因子结合和基因表达改变。在GATA2区域，PU.1，RUNX1和GATA2结合上游增强子，维持染色质环路并促进转录(图11B)，在ATRA处理后，引起转录因子结合的丧失，破坏染色质环从而抑制Gata2的转录。在Zbtb16区域中，5'和3'环维持对照细胞中Zbtb16的连续转录。在ATRA处理后，与Gata2相反，PU.1的结合导致染色质环的破坏从而抑制了Zbtb16的转录(图11C)。

Claims

1.一种样品状态转变效应基因的鉴定方法，所述效应基因的表达受到样品状态转变中染色质相互作用变化的影响，其包括下列步骤：

(2)将步骤(1)获得的差异信息建立关联，获得与状态转变中转录调控区域的染色质相互作用差异有关的基因可识别行为差异，从而鉴定所述效应基因；

其中通过以下步骤获得步骤(1)中存在于基因转录调控区域的染色质相互作用差异：

(a)使用能够标记活性启动子或增强子的组蛋白相关抗体，采用ChIP-seq技术鉴定样品基因组中处于激活状态的启动子和/或增强子的位点；

(b)鉴定所有染色质相互作用的发生区域；

(c)整合步骤(a)和步骤(b)所获得的信息，得到位于激活状态的启动子和增强子之间存在的染色质相互作用，即存在于基因转录调控区域的染色质相互作用；和

(d)将不同样品之间存在于基因转录调控区域的染色质相互作用进行比较，得到存在于基因转录调控区域的染色质相互作用差异；

其中步骤(2)通过以下步骤进行：

a)将基因表达量差异与基因调控性染色质相互作用差异结合，选择在不同状态的样品中，存在于增强或减弱的基因调控性染色质相互作用内部、同时表达量也具有显著变化的基因作为效应基因；或

b)将基因组转录调控区域中转录因子结合模体分布的差异与基因调控性染色质相互作用差异信息相结合，选择在不同状态的样品中，存在于增强的或减弱的基因调控性染色质相互作用内部、转录因子结合模体分布也发生显著变化的基因作为效应基因；或

c)将基因表达量差异、转录因子结合模体分布的差异与基因调控性染色质相互作用差异信息相结合，选择在不同状态的样品中，存在于增强的或减弱的基因调控性染色质相互作用内部、基因组转录调控区域中结合模体分布发生显著变化、同时表达量也具有显著变化的基因作为效应基因。

2.根据权利要求1所述的方法，其中所述样品是细胞。

3.根据权利要求1所述的方法，其中所述基因可识别行为差异包括基因表达量的差异和/或基因调控区域基因组序列中结合模体分布的差异。

4.根据权利要求3所述的方法，其中所述基因表达量的差异是mRNA表达量差异或蛋白质表达量差异。

5.根据权利要求1所述的方法，其中通过以下步骤获得步骤(1)中的基因可识别行为差异：

i)获得处于第一状态和第二状态的样品；

ii)取部分处于第一状态和第二状态的样品，分别进行转录表达分析，并比较样品间的mRNA表达量差异。

6.根据权利要求5所述的方法，其中所述转录表达分析采用RNA测序法。

7.根据权利要求5所述的方法，其中进一步包括：

iii)取部分处于第一状态和第二状态的样品，分别进行染色质开放区域序列的分析，并分析染色质开放区域序列中所分布的转录因子结合模体，进一步比较处于第一状态和第二状态的样品间的转录因子结合模体的分布差异。

8.根据权利要求7所述的方法，其中通过以下步骤分析步骤(1)中存在于基因转录调控区域的染色质相互作用差异：

iv)取部分处于第一状态和第二状态的样品，鉴定分别处于激活状态的启动子和/或增强子信息；

v)另取部分处于第一状态和第二状态的样品，采用高通量染色质构象捕获技术，获得全基因组染色质相互作用的信息；或者利用4C或5C法获得局部染色质相互作用的信息；

vi)将参考基因组序列划分成一定大小的区域，所述区域大小在1-40kb之间，基于步骤iv)获得的活性状态的启动子和增强子信息，通过比对分别获得包含有活性的启动子和增强子位点的区域，所述包含启动子的区域命名为基因区域，将包含增强子序列的区域命名为调控区域；

随后，结合步骤v)获得的染色质相互作用的频率信号，识别出发生于基因区域和调节区域之间的染色质相互作用频率信号，从而得到基因调控相关染色质相互作用；然后将处于第一状态和第二状态的样品之间的基因调控相关染色质相互作用进行比较，其中具有统计显著性差异的，被鉴定为基因调控性染色质相互作用差异(differential gene-regulatory interaction)，包括相对于处于第一状态的样品，处于第二状态的样品中增强的基因调控性染色体相互作用和/或减弱的基因调控性染色体相互作用。

9.根据权利要求8所述的方法，其中所述步骤iv)中的鉴定采用ChIP-seq法，所述ChIP-seq法中所使用的抗体为H3K4me3和H3K27ac的结合抗体，所述抗体分别结合H3K4me3和H3K27ac形成信号峰，分别代表了处于激活状态的启动子和增强子位点。

10.根据权利要求8所述的方法，其中所述高通量染色质构象捕获技术是Hi-C法、insitu Hi-C法、BL-Hi-C法或ChIA-PET法。

11.根据权利要求1所述的方法，此外还包括对筛选出来的效应基因进行功能研究以确定其功能的步骤。

12.根据权利要求1所述的方法，其中所述样品状态转变通过下述方式实现：化学试剂诱导、自然分化和/或物理刺激。