基于ChIP-seq数据HMM方法识别全基因组的差异组蛋白修饰位点

1、基于ChIP-seq数据HMM方法识别全基因组的差异组蛋白修饰位点摘要目的：表观遗传修饰是调控基因表达和基因组功能的一个主要因素。在不同的表观遗传修饰中，差异组蛋白修饰位点（DHMSs）是不同细胞类型、时期和环境影响时，表观遗传动态性质和基因表达调控的一个研究热点。为了测定全基因组的组蛋白修饰，ChIP-seq技术是一种有效的方法。因此，通过比较两个ChIP-seq文库可以识别潜在的DHMSs。结果：我们的目的是识别DHMSs，提出一种称为ChIPDiff的方法来通过ChIP-seq测定的数据全基因组比对组蛋白修饰位点。基于观察的ChIP片段数，提出了一个隐马模型的方法推断每个基因组位置的组蛋

2、白修饰变化状态。我们通过比对小鼠ESC和NPC细胞的H3K27me3修饰位点来评估ChIPDiff的效果。我们证明了此方法确定H3K27me3 的DHMSs具有高灵敏度，特异性和重复性。进一步应用ChIPDiff揭示不同细胞时期的差异H3K4me3和H3K36me3位点。我们研究中的比对有很多有趣的生物学发现。1.介绍真核DNA是被打包到一个由周围环绕组蛋白的DNA的重复核小体组成的染色质结构。组蛋白可以发生大量的翻译后修饰如，甲基化，乙酰化，磷酸化和泛素化。组蛋白修饰影响基因表达和基因组功能。大量实验证明一些组蛋白甲基化类型在生物学过程中起主要作用。一个典型的例子是在哺乳动物胚胎干细胞通过H

3、3K27me3抑制发育调控维持干细胞多能性。在癌症中也特异的发现一些表观遗传K27干细胞标记。此外，H3K9me3、H3K9me2和癌细胞中沉默肿瘤抑制基因相关。因此，特异基因组位置的差异组蛋白修饰密度，文中称为差异组蛋白修饰位点“DHMS”，在不同细胞类型，时期和环境影响是比较研究的重点。我们可以用染色质免疫共沉淀（ChIP）来测定组蛋白修饰信号，抗体用于富集修饰位点的DNA片段。在过去的几年开发了几种基于ChIP的技术，包括ChIP-chip, ChIP-PET and ChIP-SAGE，用于大规模基因组区域的组蛋白修饰和转录因子结合位点研究。随着最近超高通量测序技术如Illumina/

4、 Solexa GA 测序的产生，ChIP-seq成为一个主要的高覆盖、高分辨率和低成本的方法。ChIP-seq的基本思想是读取ChIP富集的序列的一端，接着映射这些短读称为tag到基因组上以找到这些片段的基因组位置。一个ChIP文库中有百万个tag标签测序，形成一个代表全基因组与组蛋白修饰位点和转录因子结合位点的ChIP片段数的谱。受到ChIP-seq在单个文库识别组蛋白修饰位点的鼓舞，我们想是否可以通过计算的比较不同细胞类型和实验条件的两条ChIP-seq文库来识别DHMS。Mikkelsen等人测定了小鼠ESC、NPC和MEF细胞的H3K4me3 (K4) 和 K27位点，比较三种类型启

5、动子区域修饰位点的发生。他们研究的局限在于修饰位点是定性的比较而非定量。一个例子说明了这种局限，K4调控K1f4,已知其和基因表达正相关。K1f4在ESC和NPC启动子定性分析中都标记K4，因此不能解释在ESC的K1f4上调。另一方面，定量比较表明ESC的K1f4启动子的K4密度比NPC多5倍，这和表达变化是一致的。据我们所知，几乎没有全基因组定量比较两个ChIP-seq文库的文献。受芯片分析的启发，一个简单的解决这个问题的方法是将基因组分为箱bins，计算每个binChIP片段数的倍数变化。然而，fold-change方法对由ChIP片段随机样本的技术变化时敏感的。本文中，我们提出的方法称为

6、ChIPDiff通过考虑连续bin之间的相关性改进了fold-change方法。我们用隐马模型建立相关性，转移概率用一种无监督方式自动训练。接下来通过训练HMM参数来推断组蛋白修饰状态的变化。为了评估ChIPDiff的性能，我们首先比较Mikkelsen数据ESC和NPC的K27文库。在全基因组识别了4277个k27的DHMS区域。三个标准显示效果是令人满意的：（a）敏感性：2006年在高度保守的非编码元件中，80%的从基因表达推断的DHMSs被ChIPDiff确定。（b）特异性：基于非细胞特异性控制比对，我们估计识别的DHMS区域的假阳性率是0.19%。（c）重复度：检查两个独立的子集的结果

7、的交集，显示3-4百万个tags测序的57.4的DHMSs在技术上重现，评价结果还表明，在所有三个方面的定性分析，该方法优于fold-change的方法。我们进一步应用ChIPDiff到H3K4me3（K4）和H3K36me3（K36），发现这两种类型组蛋白修饰的DHMSs和研究了他们在干细胞分化潜在的生物的作用。研究中有几个有趣的生物学发现。2.方法2.1确定组蛋白修饰位点给定来个ChIP-seq文库，L1和L2，识别DHMSs的第一步是确定L1和L2的组蛋白修饰假定的位点。这部分详述这一步。ChIP-seq实验产生的原始数据的tags被映射到基因组，获得它们的位置和方向。由于ChIP-se

8、q实验的PCR过程，大量的tags可能源于一个单一的ChIP片段。为了移除这一重复性，映射到相同位置和相同方向的tags被作为一个单一的copy。注意到在ChIP-seq协议一个单一的tag是通过测序一个ChIP片段的末端得到的，平均长度是200bp。因此我们通过其方向的100bp转移tag的位置近似估计响应ChIP片段的中心。全基因组被分成1k-bp的bin，计算每个bin的ChIP片段中心数。预处理过程之后，产生ChIP片段数谱。考虑到基因组有m个bin，谱L1和L2分别表示为X1 = x1.1，x1.2，.x1.m 和X2 = x2.1，x2.2，.x2.m。其中xij是在Li中第j个b

9、in的片段数。为了描述每个bin中片段的结合富集，我们定义F值标准化测序的深度：其中n1和n2是L1和L2测序片段的总数，如图。 Mikkelson et al.(2007)和Robertson et al.(2007)指出有与重复序列区域的存在，并不是所有的bin都能在tag映射程序中检测到。让记为基因组“有效”的bin，分值F的期望在有效bin时是F(i)/(m×)，等于2/(m×)。Mikkelson et al.(2007)估计小鼠基因组的等于0.7。如果一个bin的F值大于2/(m×)，我们标记其为一个推测的组蛋白修饰位点。1k bp内的连续修饰位点彼此

10、分开被合并为组蛋白修饰区域。2.2用Fold-change方法定量的比较修饰强度为了便于定义和描述，文章其他部分将介绍的基于推定的组蛋白修饰区域在2.1介绍，假设一个区域包含k个bin，我们定义L1和L2的ChIP片段数分别为x1.i，x2.i，在区域的第i个bin（i=1，1，k）。组蛋白修饰表现出对各种动力性和化学计量性。对一个ChIP实验，我们定义文库Lj的第i个bin的修饰强度是任意ChIP片段来自ChIP过程第i个bin的概率，定义为pj，i。由于提取和测序ChIP片段是一个随机抽样过程，文库Lj的第i个bin的观察片段xj，i的后验概率，强度的条件概率pj，i，近似服从二项分布：（

11、1）我们接下来估计先验概率pj，i服从beta分布：（2）B(,)是beta函数。注意到beta分布先于二项是共轭的，所以条件概率也服从beta分布，期望等于。在我们的应用中，参数和设为1和m，m是基因组中bin的总数（详见补充方法）。我们定义一个DHMS，当一个bin内L1和L2的强度比值大于（L1富集DHMS）或者小于1/(L2富集DHMS)。是一个预先确定的阈值，值1。一个简单识别DHMSs的方法是估计ChIP片段数的期望强度（更好的是对数比）的倍数变化，如下：（3）基于方程（3）的对数比估计显示图1（a）。fold-change法的一个缺陷是由于随机抽样引起技术差异。图1（b）显示一个

12、RI-plot描述了依据强度的log比值变化。当强度相对较小，log值的变化太高，这可能引起大量的假阳性。2.3一个基于隐马模型的方法识别DHMSs组蛋白修饰通常发生在连续区域范围是几百甚至上千个核苷酸。因此可以期望连续的bin测量的强度变化可能强相关。通过观察ChIP-seq谱支持这一观点。例如，图1（a）的log比值谱的自相关是0.84。在ChIP-chip数据分析中，Li et al.(2005)年设计的HMM模型构建连续探针之间的信号相关成功的应用于识别p53结合位点，表示HMM在我们研究中应用的潜在可能性。在此我们提出一个基于HMM的方法，ChiPDiff来解决这一问题。我们定义Si

13、为第i个bin的组蛋白修饰变化状态（i=1到k），基于2.2对于DHMS的定义，状态Si为以下三个值之一：0：无差别位点，if 1/p1，i/p2，i ；1：L1富集DHMS，if p1，i/p2，i；2：L2富集DHMS，if p1，i/p2，i1/。我们建模bin间的相关性作为一个一阶马尔可夫链Pr( Si|S0,S1，., Si-1)= Pr(Si|Si-1)，S0是区域内第一个bin前的起始状态。一个HMM实施是通过观察片段数推断状态的后验概率分布。HMM的三个特征:起始状态S0的先验概率，emission发射概率，和状态转移概率。初始状态S0采用固定值0，因为我们假定两个文库中区域起

14、始位置是组蛋白修饰缺乏的基因组位置。我们通过整合所有可能的Si值的p1，i和p2，i得到emission发射概率读者可以参考补充方法的详细推导。在等式（4）中，服从二项分布（1），服从beta分布（2）。转移概率列表由Baum-Welch算法训练得到，采用期望最大化（EM）步骤以无监督的方式从隐藏状态迭代估计HMM的参数。训练过程中，传输参数初始化是统一的，初始状态S0和状态传输概率如以上描述确定。因为转移概率表在整个基因组是相同的，是通过所有推定的组蛋白修饰区域转移频率累加训练的（train）。在ChiPDiff的最后一步，每个bin中的概率分布状态由forward-backward算法推断

15、。如果bin的后验概率大于置信阈值（0<<1)当Si=1或Si=2定为一个DHMS区。连续的没有缝隙的DHMS被合并为一个DHMS。ChiPDiff最大计算量的一步是训练转移概率表。两个策略可以减少计算量(a) 训练HMM之前，发射概率的积分被数值计算的而且被编写成一张查询列表。(b)我们允许转移概率列表基于从推定组蛋白修饰区域随机选择子集训练。3.结果我们应用ChIPDiff处理Mikkelson实验数据，ChIPDiff的的性能通过比较小鼠ESC和NPC的H3K27me3文库评估。我们又应用ChIPDiff处理H3K4me3和H3K36me3数据发现了DHMSs而且研究它们在干

16、细胞分化中潜在的生物学作用。3.1 H3K27me3数据评估选用H3K27me3评估的原因是因为它的DHMSs在高度保守的非编码元件（HCNEs）已经有人研究。而且，K27优先标记基因区域功能作为抑制子，这有利于我们利用表达数据间接的验证。我们用ChIPDiff比较ESC和NPC的K27文库，fold-change阈值设为3.0置信阈值为0.95.HMM随机训练10000次选定组蛋白修饰区域，26230bins认定为DHMS是，对应于4722连续区域。它们中3,833 (81.2%)区域ESC富集，889 (18.8%)NPC富集，这意味着细胞分化时期K27消耗的整体趋势。我们首次评估了ChI

17、PDiff的性能通过确定其生物学意义，如敏感性。Bernstein发现K27在ESC中富集在高度保守的非编码元件（HCNEs），抑制发育调控子的数量来维持细胞的stemness。这些组蛋白标记在不同分化细胞中消失。HCNEs中，我们选择了223个基因，Mikkelson研究了它们的表达。因为K27作为功能抑制子，这些中的一些被K27标记的HCNEs在NPC中上调，我们认为在这些基因DHMSs被确定。与预期相同，一个包含30个上调基因的子集被确定，标准化超过4倍。它们中80%被标记的由ChIPDiff识别的DHMSs在启动子区域转录起始位点±1kb。相反，193个基因中只有19.2%在

18、NPC中被DHMSs标记的不上调。为了检验ChIPDiff的特异性，我们需要评估错误识别的不是细胞特异的DHMS区域的片段。针对这一目的，我们将这一文库分为两个技术复制本：Lesc，rep1和Lesc，rep2，Lnpc，rep1和Lnpc，rep2。复制样本的tag组成取自ChIP-seq实验的不同通道，有相似的测序depth（见补充表四复制本的通道）。通过合并tags产生两个新的文库分别是Lesc，rep1和Lnpc，rep1，Lesc，rep2和Lnpc，rep2。因为复制本有相似的测序depth，两个文库的差异不是细胞特异的可能只是实验技术变化的影响。比较这些非细胞特异的控制集，Chi

19、PDiff识别出9个差异的区域，因此我们估计在识别细胞特异比较时DHMS区时假阳性率为0.19（9/4722）。我们通过构建两个独立的细胞特异比对途径passes检测重复性：Lesc，rep1vs.Lnpc，rep1，Lesc，rep2vs.Lnpc，rep2。为了测量重复性，我们定义一个分数作为两个passes均识别的DHMSs数目与单个pass识别DHMSs平均值的比值。结果得ChIPDiff重复分数为57.4%。注意到重复性的条件是重复本的测序depth在3到4百万个tags（补充表4）。为了比较不同方法的效果，我们重复比较fold-change法和定量方法的敏感性、特异性和重复度。在定

20、量方法中，ESC和NPC的K27修饰位点单个识别用Mikkelson的方法，K27位点只在单个细胞类型标记识别为DHMSs。连续的DHMSs合并为DHMS区。为了公正的比对，阈值调整使所有三种方法DHMSs区域相似数目被确定（因为阈值选的离散值，所以这个数目不同）。评估结果总结见表1。ChIPDiff所有三方面均优于其他两种方法。Fold-change和定量的方法都有高的假阳性率，表示这些方法对技术变化和实验偏差是敏感的。3.2应用H3K4me3和H3K36me3数据我们扩展我们的数据研究到H3K4me3和H3K36me3。这两种修饰类型以不同的方式正向调控基因表达。Guenther2007发

21、现K4在基因转录起始标记活性启动子，而K36发生在基因区作为延伸的标志。我们之前的研究（Zhao等人）也显示K4和K27在人类ESC中建立不同的基因组区域活性和非活性染色质结构。这引起我们研究NPC和ESC之间这些组蛋白修饰的DHMSs的兴趣。此外，K4位点通常在ChIP-seq谱转录起始位点周围急剧以点状的模式出现，而K36位点出现以更宽的模式，提供一种全面的检验-bed评估我们方法对于不同组蛋白修饰的适用性。我们用3.1提到的ChiPDiff运行相同的文库。结果见表2。连续的DHMSs合并在一起。值得注意的是，K4在ESC富集的DHMSs远远大于在K4在NPC富集的DHMSs。考虑到这种不

22、平衡也出现在K27，我们假设被K4和K27标记的二价染色质结构可能相关（Bernstein2006）。接下来的分析中，我们发现3833个ESC的K27 DHMSs中的1961（51.2%）个与K4的DHMSs重叠。相反，K36和K27倾向于互相排斥：只有8个（0.21%）个DHMS重叠。为了研究DHMSs和基因表达的相关性，我们注释DHMSs区的Refseq基因和基因表达数据。为了去除冗余，基因注释最长ORF选择是如果相同的转录本注释到相同的基因，结果筛选一共18795个唯一的基因。如图2所示，K4和K36共调控基因表达强显著相关。这与之前的研究相一致。1085个基因在ESC中上调，791（72.9%）和ESC富集K4和K36相关，表示基因表达可能潜在的由DHMS预测。显然，ESC两个关键的转录因子，Nanog和Oct4，由K4和K36的DHMSs标记，也暗示这些ESC的组蛋白修饰标记在影响转录调控网络其主