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

《基于ChIP-seq数据HMM方法识别全基因组的差异组蛋白修饰位点》由会员分享，可在线阅读，更多相关《基于ChIP-seq数据HMM方法识别全基因组的差异组蛋白修饰位点（6页珍藏版）》请在金锄头文库上搜索。

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

2、的数据全基因组比对组蛋白修饰位点。基于观察的 ChIP 片段数，提出 了一个隐马模型的方法推断每个基因组位置的组蛋白修饰变化状态。我们通过 比对小鼠 ESC 和 NPC 细胞的 H3K27me3 修饰位点来评估 ChIPDiff 的效果。我们 证明了此方法确定 H3K27me3 的 DHMSs 具有高灵敏度，特异性和重复性。进一 步应用 ChIPDiff 揭示不同细胞时期的差异 H3K4me3 和 H3K36me3 位点。我们研 究中的比对有很多有趣的生物学发现。1.1.介绍介绍真核 DNA 是被打包到一个由周围环绕组蛋白的 DNA 的重复核小体组成的染色质 结构。组蛋白可以发生大量的翻译后修

3、饰如，甲基化，乙酰化，磷酸化和泛素 化。组蛋白修饰影响基因表达和基因组功能。大量实验证明一些组蛋白甲基化 类型在生物学过程中起主要作用。一个典型的例子是在哺乳动物胚胎干细胞通 过 H3K27me3 抑制发育调控维持干细胞多能性。在癌症中也特异的发现一些表观 遗传 K27 干细胞标记。此外，H3K9me3、H3K9me2 和癌细胞中沉默肿瘤抑制基因 相关。因此，特异基因组位置的差异组蛋白修饰密度，文中称为差异组蛋白修 饰位点“DHMS”，在不同细胞类型，时期和环境影响是比较研究的重点。我们 可以用染色质免疫共沉淀（ChIP）来测定组蛋白修饰信号，抗体用于富集修饰 位点的 DNA 片段。在过去的几

4、年开发了几种基于 ChIP 的技术，包括 ChIP-chip, ChIP-PET and ChIP-SAGE，用于大规模基因组区域的组蛋白修饰和转录因子结 合位点研究。随着最近超高通量测序技术如 Illumina/ Solexa GA 测序的产生， ChIP-seq 成为一个主要的高覆盖、高分辨率和低成本的方法。ChIP-seq 的基本 思想是读取 ChIP 富集的序列的一端，接着映射这些短读称为 tag 到基因组上以 找到这些片段的基因组位置。一个 ChIP 文库中有百万个 tag 标签测序，形成一 个代表全基因组与组蛋白修饰位点和转录因子结合位点的 ChIP 片段数的谱。受 到 ChIP-

5、seq 在单个文库识别组蛋白修饰位点的鼓舞，我们想是否可以通过计算 的比较不同细胞类型和实验条件的两条 ChIP-seq 文库来识别 DHMS。Mikkelsen 等人测定了小鼠 ESC、NPC 和 MEF 细胞的 H3K4me3 (K4) 和 K27 位点，比较三种 类型启动子区域修饰位点的发生。他们研究的局限在于修饰位点是定性的比较 而非定量。一个例子说明了这种局限，K4 调控 K1f4,已知其和基因表达正相关。 K1f4 在 ESC 和 NPC 启动子定性分析中都标记 K4，因此不能解释在 ESC 的 K1f4 上调。另一方面，定量比较表明 ESC 的 K1f4 启动子的 K4 密度比

6、NPC 多 5 倍， 这和表达变化是一致的。 据我们所知，几乎没有全基因组定量比较两个 ChIP-seq 文库的文献。受芯片分 析的启发，一个简单的解决这个问题的方法是将基因组分为箱 bins，计算每个binChIP 片段数的倍数变化。然而，fold-change 方法对由 ChIP 片段随机样本 的技术变化时敏感的。本文中，我们提出的方法称为 ChIPDiff 通过考虑连续 bin 之间的相关性改进了 fold-change 方法。我们用隐马模型建立相关性，转 移概率用一种无监督方式自动训练。接下来通过训练 HMM 参数来推断组蛋白修 饰状态的变化。 为了评估 ChIPDiff 的性能，我们

7、首先比较 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 的方法。

8、 我们进一步应用 ChIPDiff 到 H3K4me3（K4）和 H3K36me3（K36），发现这两种 类型组蛋白修饰的 DHMSs 和研究了他们在干细胞分化潜在的生物的作用。研究 中有几个有趣的生物学发现。2.2.方法方法2.1 确定组蛋白修饰位点 给定来个 ChIP-seq 文库，L1 和 L2，识别 DHMSs 的第一步是确定 L1 和 L2 的组 蛋白修饰假定的位点。这部分详述这一步。 ChIP-seq 实验产生的原始数据的 tags 被映射到基因组，获得它们的位置和方 向。由于 ChIP-seq 实验的 PCR 过程，大量的 tags 可能源于一个单一的 ChIP 片 段。为了移除

9、这一重复性，映射到相同位置和相同方向的 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

10、 中第 j 个 bin 的片段数。为了描述每个 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)，我们标 记其为

11、一个推测的组蛋白修饰位点。1k bp 内的连续修饰位点彼此分开被合并 为组蛋白修饰区域。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。由

12、于提取和测序 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)。 是一个预先确定的阈值，

13、值1。一个 简单识别 DHMSs 的方法是估计 ChIP 片段数的期望强度（更好的是对数比）的倍数变化，如下：（3）基于方程（3）的对数比估计显示图 1（a）。 fold-change 法的一个缺陷是由于随机抽样引起技术差异。图 1（b）显示一个 RI-plot 描述了依据强度的 log 比值变化。当强度相对较小，log 值的变化太高， 这可能引起大量的假阳性。 2.3 一个基于隐马模型的方法识别 DHMSs 组蛋白修饰通常发生在连续区域范围是几百甚至上千个核苷酸。因此可以期望 连续的 bin 测量的强度变化可能强相关。通过观察 ChIP-seq 谱支持这一观点。 例如，图 1（a）的 log

14、 比值谱的自相关是 0.84。在 ChIP-chip 数据分析中，Li et al.(2005)年设计的 HMM 模型构建连续探针之间的信号相关成功的应用于识 别 p53 结合位点，表示 HMM 在我们研究中应用的潜在可能性。在此我们提出一 个基于 HMM 的方法，ChiPDiff 来解决这一问题。 我们定义 Si 为第 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

15、，i1/。 我们建模 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），服从 b

16、eta 分布（2）。转移概率列表由 Baum-Welch 算法训练得到，采用期望最大化（EM）步骤以无监 督的方式从隐藏状态迭代估计 HMM 的参数。训练过程中，传输参数初始化是统 一的，初始状态 S0 和状态传输概率如以上描述确定。因为转移概率表在整个基 因组是相同的，是通过所有推定的组蛋白修饰区域转移频率累加训练的 （train）。 在 ChiPDiff 的最后一步，每个 bin 中的概率分布状态由 forward-backward 算 法推断。如果 bin 的后验概率大于置信阈值 （01)当 Si=1 或 Si=2 定为一个 DHMS 区。连续的没有缝隙的 DHMS 被合并为一个 DHMS。ChiPDiff 最大 计算量的一步是训练转移概率表。两个策略可以减少计算量(a) 训练 HMM 之前， 发射概率的积分被数值计算的而且被编写成一张查询列表。(b)我们允许转移概 率列表基于从推定组蛋白修饰区域随机选择子集训练。3.3.结果结果我们应用 ChIPDiff 处理 Mikkel