搜索
基于ChIP-seq数据HMM方法识别全基因组的差异组蛋白修饰位点摘要目的：表观遗传修饰是调控基因表达和基因组功能的一个主要因素。在不同的表观遗传修饰中，差异组蛋白修饰位点（DHMSs）是不同细胞类型、时期和环境影响时，表观遗传动态性质和基因表达调控的一个研究热点。为了测定全基因组的组蛋白修饰，ChIP-seq技术是一种有效的方法。因此，通过比较两个ChIP-seq文库可以识别潜在的DHMSs。结果：我们的目的是识别DHMSs，提出一种称为ChIPDiff的方法来通过ChIP-seq测定的数据全基因组比对组蛋白修饰位点。基于观察的ChIP片段数，提出了一个隐马模型的方法推断每个基因组位置的组蛋白修饰变化状态。我们通过比对小鼠ESC和NPC细胞的H3K27me3修饰位点来评估ChIPDiff的效果。我们证明了此方法确定H3K27me3的DHMSs具有高灵敏度，特异性和重复性。进一步应用ChIPDiff揭示不同细胞时期的差异H3K4me3和H3K36me3位点。我们研究中的比对有很多有趣的生物学发现。1.介绍真核DNA是被打包到一个由周围环绕组蛋白的DNA的重复核小体组成的染色质结构。组蛋白可以发生大量的翻译后修饰如，甲基化，乙酰化，磷酸化和泛素化。组蛋白修饰影响基因表达和基因组功能。大量实验证明一些组蛋白甲基化类型在生物学过程中起主要作用。一个典型的例子是在哺乳动物胚胎干细胞通过H3K27me3抑制发育调控维持干细胞多能性。在癌症中也特异的发现一些表观遗传K27干细胞标记。此外，H3K9me3、H3K9me2和癌细胞中沉默肿瘤抑制基因相关。因此，特异基因组位置的差异组蛋白修饰密度，文中称为差异组蛋白修饰位点“DHMS”，在不同细胞类型，时期和环境影响是比较研究的重点。我们可以用染色质免疫共沉淀（ChIP）来测定组蛋白修饰信号，抗体用于富集修饰位点的DNA片段。在过去的几年开发了几种基于ChIP的技术，包括ChIP-chip,ChIP-PETandChIP-SAGE，用于大规模基因组区域的组蛋白修饰和转录因子结合位点研究。随着最近超高通量测序技术如Illumina/SolexaGA测序的产生，ChIP-seq成为一个主要的高覆盖、高分辨率和低成本的方法。ChIP-seq的基本思想是读取ChIP富集的序列的一端，接着映射这些短读称为tag到基因组上以找到这些片段的基因组位置。一个ChIP文库中有百万个tag标签测序，形成一个代表全基因组与组蛋白修饰位点和转录因子结合位点的ChIP片段数的谱。受到ChIP-seq在单个文库识别组蛋白修饰位点的鼓舞，我们想是否可以通过计算的比较不同细胞类型和实验条件的两条ChIP-seq文库来识别DHMS。Mikkelsen等人测定了小鼠ESC、NPC和MEF细胞的H3K4me3(K4)和K27位点，比较三种类型启动子区域修饰位点的发生。他们研究的局限在于修饰位点是定性的比较而非定量。一个例子说明了这种局限，K4调控K1f4,已知其和基因表达正相关。K1f4在ESC和NPC启动子定性分析中都标记K4，因此不能解释在ESC的K1f4上调。另一方面，定量比较表明ESC的K1f4启动子的K4密度比NPC多5倍，这和表达变化是一致的。据我们所知，几乎没有全基因组定量比较两个ChIP-seq文库的文献。受芯片分析的启发，一个简单的解决这个问题的方法是将基因组分为箱bins，计算每个binChIP片段数的倍数变化。然而，fold-change方法对由ChIP片段随机样本的技术变化时敏感的。本文中，我们提出的方法称为ChIPDiff通过考虑连续bin之间的相关性改进了fold-change方法。我们用隐马模型建立相关性，转移概率用一种无监督方式自动训练。接下来通过训练HMM参数来推断组蛋白修饰状态的变化。为了评估ChIPDiff的性能，我们首先比较Mikkelsen数据ESC和NPC的K27文库。在全基因组识别了4277个k27的DHMS区域。三个标准显示效果是令人满意的：（a）敏感性：2006年在高度保守的非编码元件中，80%的从基因表达推断的DHMSs被ChIPDiff确定。（b）特异性：基于非细胞特异性控制比对，我们估计识别的DHMS区域的假阳性率是0.19%。（c）重复度：检查两个独立的子集的结果的交集，显示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-seq实验的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个bin的片段数。为了描述每个bin中片段的结合富集，我们定义F值标准化测序的深度：其中n1和n2是L1和L2测序片段的总数，如图。Mikkelsonetal.(2007)和Robertsonetal.(2007)指出有与重复序列区域的存在，并不是所有的bin都能在tag映射程序中检测到。让η记为基因组“有效”的bin，分值F的期望在有效bin时是ΣF(i)/(m×η)，等于2/(m×η)。Mikkelsonetal.(2007)估计小鼠基因组的η等于0.7。如果一个bin的F值大于2/(m×η)，我们标记其为一个推测的组蛋白修饰位点。1kbp内的连续修饰位点彼此分开被合并为组蛋白修饰区域。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，近似服从二项分布：（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）显示一个RI-plot描述了依据强度的log比值变化。当强度相对较小，log值的变化太高，这可能引起大量的假阳性。2.3一个基于隐马模型的方法识别DHMSs组蛋白修饰通常发生在连续区域范围是几百甚至上千个核苷酸。因此可以期望连续的bin测量的强度变化可能强相关。通过观察ChIP-seq谱支持这一观点。例如，图1（a）的log比值谱的自相关是0.84。在ChIP-chip数据分析中，Lietal.(2005)年设计的HMM模型构建连续探针之间的信号相关成功的应用于识别p53结合位点，表示HMM在我们研究中应用的潜在可能性。在此我们提出一个基于HMM的方法，ChiPDiff来解决这一问题。我们定义Si为第i个bin的组蛋白修饰变化状态（i=1到k），基于2.2对于DHMS的定义，状态Si为以下三个值之一：α0：无差别位点，if1/τ≤p1，i/p2，i≤τ；α1：L1富集DHMS，ifp1，i/p2，i＞τ；α2：L2富集DHMS，ifp1，i/p2，i＜1/τ。我们建模bin间的相关性作为一个一阶马尔可夫链Pr(Si|S0,S1，...,Si-1)=Pr(Si|Si-1)，S0是区域内第一个bin前的起始状态。一个HMM实施是通过观察片段数推断状态的后验概率分布。HMM的三个特征:起始状态S0的先验概率，emission发射概率，和状态转移概率。初始状态S0采用固定值α0，因为我们假定两个文库中区域起始位置是组蛋白修饰缺乏的基因组位置。我们通过整合所有可能的Si值的p1，i和p2，i得到emission发射概率读者可以参考补充方法的详细推导。在等式（4）中，服从二项分布（1），服从beta分布（2）。转移概率列表由Baum-Welch算法训练得到，采用期望最大化（EM）步骤以无监督的方式从隐藏状态迭代估计HMM的参数。训练过程中，传输参数初始化是统一的，初始状态S0和状态传输概率如以上描述确定。因为转移概率表在整个基因组是相同的，是通过所有推定的组蛋白修饰区域转移频率累加训练的（train）。在ChiPDiff的最后一步，每个bin中的概率分布状态由forward-backward算法推断。如果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而且研究它们在干细胞分化中潜在的生物学作用。3.1H3K27me3数据评估选用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消耗的整体趋势。我们首次评估了ChIPDiff的性能通过确定其生物学意义，如敏感性。Bernstein发现K27在ESC中富集在高度保守的非编码元件（HCNEs），抑制发育调控子的数量来维持细胞的stemness。这些组蛋白标记在不同分化细胞中消失。HCNEs中，我们选择了223个基因，Mikkelson研究了它们的表达。因为K27作为功能抑制子，这些中的一些被K27标记的HCNEs在NPC中上调，我们认为在这些基因DHMSs被确定。与预期相同，一个包含3