我们的贝叶斯变点（BCP）算法基于利用隐马尔可夫模型对变点建模的前沿进展，并将其应用于染色质免疫沉淀测序（ChIPseq）数据分析。BCP在宽峰和尖峰数据类型中均表现良好，尤其擅长准确识别弥散型组蛋白富集的稳定且可重复的区域。
需要JoVE订阅才能观看此内容。 请登录或开始免费试用
本文内容
摘要
摘要
ChIPseq 是一种广泛用于研究蛋白质-DNA 相互作用的技术。通过下一代测序对蛋白质结合的 DNA 进行测序，并将短序列读段比对到参考基因组，从而生成读段密度谱。富集区域表现为峰，其峰形通常因靶蛋白的不同而存在显著差异1。例如，转录因子通常以位点和序列特异性的方式结合，往往产生尖锐的峰；而组蛋白修饰则更为广泛，其特征是出现宽广、弥散的富集岛2。可靠地识别这些区域正是我们工作的重点。
用于分析ChIPseq数据的算法采用了多种方法，从启发式方法到3-5 到更严格的统计模型， 例如 隐马尔可夫模型（HMMs）6-8我们寻求一种解决方案，以最大限度地减少对难以定义的临时参数的依赖，这些参数往往会降低分辨率，并削弱工具的直观可用性。针对基于隐马尔可夫模型（HMM）的方法，我们的目标是简化参数估计流程，并避免常被采用的简单有限状态分类。
此外，传统的ChIPseq数据分析需要将预期的读段密度分布特征归类为点状或弥散状，然后分别应用相应的分析工具。我们进一步旨在用一种更通用的单一模型取代上述两种独立模型，该模型能够有效应对各种类型的数据。
为实现这些目标，我们首先构建了一个统计框架，该框架利用隐马尔可夫模型（HMM）领域的前沿进展9，以自然的方式对ChIP-seq数据结构进行建模，且仅使用显式公式——这一创新对其性能优势至关重要。相较于启发式模型，我们的HMM通过贝叶斯模型容纳了无限隐状态。我们将该模型应用于识别读段密度中的合理变点，从而进一步定义富集区域的片段。分析结果表明，我们的贝叶斯变点（BCP）算法具有较低的计算复杂度，表现为运行时间缩短和内存占用减少。BCP算法在尖峰状峰（punctate peak）和弥散岛状区域（diffuse island）识别中均成功应用，具备稳健的准确性且所需用户自定义参数较少。这体现了该算法的多功能性和易用性。因此，我们认为该方法可便捷地应用于多种数据类型和广泛的终端用户，且结果易于进行比较与对照，使其成为一种有助于不同研究团队之间协作与验证的优秀ChIP-seq数据分析工具。本文中，我们通过将BCP应用于已有的转录因子10,11和表观遗传数据12，展示了其实际应用价值。
方案
1. 为 BCP 分析准备输入文件
- 使用首选的短序列比对软件，将测序运行产生的短序列读段（ChIP 和输入文库）比对至相应的参考基因组。比对到的基因组位置应转换为六列的浏览器可扩展数据（BED）格式13（UCSC基因组浏览器，http://genome.ucsc.edu/），每条比对读段对应一行以制表符分隔的数据，包含比对到的染色体、起始位置（0-based）、终止位置（半开区间）、读段名称、分值（可选）以及链方向。
2a. 弥散型读段分布：弥散数据中富集区域检测的ChIP读段密度预处理
- 将 ChIP 和输入的比对位置扩展至预设的片段长度， 即DNA 酶切或超声处理过程中目标片段大小，通常约为 200 bp。随后将相邻区段内的片段计数进行汇总。默认情况下，区段大小设为估计的片段长度 200 bp。
- 在一组读段数相同的分箱中，任何可能的变点极有可能位于最外侧的边界处。因此，在两个读段数相同的相邻分箱之间的内部边界上出现变点的可能性极低。应将读段数相同的相邻分箱合并为一个单一的区块。 即bedGraph 格式13.
2b. 点状读段分布：点状数据中峰区域检测的ChIP与Input BED文件预处理
- 分别对正链和负链的ChIP读段进行重叠读段的合并。链特异性的读段密度应形成正链和负链峰的双峰分布。选择信号最富集的正链/负链峰对，以其峰顶之间的距离作为文库片段长度的估计值。
- 将ChIP和input读段向中心方向移动半个片段长度，并重新计算移动后合并的正链与负链读段的密度。该片段长度估计方法借鉴自Zhang, et al.3。具有相同合并计数的位点应归为同一区块，类似于步骤2a.2。
3. 使用我们的 BCMIX 近似方法估计每个区块的后验均值读段密度
- 将每个区块的读取密度建模为泊松分布Pois(θt)，其均值参数服从伽马分布Γ(α,β)的混合分布，并假设在任意区块边界处出现变点的先验概率为p。将Pois(θt)以G(α,β)为条件，实质上使该模型成为无限状态隐马尔可夫模型（HMM）。使用最大后验似然法估计超参数α、β和p。
- 显式计算每个区块的贝叶斯估计值θt，即E(θt|γZ)。采用计算效率更高的有限复杂度混合近似方法，替代隐马尔可夫模型中传统但耗时的前向-后向滤波算法，以估计后验均值θc。所得后验均值将被“平滑”为近似的分段常数轮廓，因此具有相同θc的区块应进一步合并，并更新其边界坐标。
4a. 弥散型读段分布：将后验均值后处理为弥散富集区段
- 使用每个新样本的输入读段数量θc 将泊松分布 Pois(λ) 作为背景速率a）并基于ChIP后验均值进行简单的假设检验以确定富集程度 θc，超过某个阈值 δ。90th-分位数是默认的d值，在大多数情况下均适用。
- 合并相邻 θc 超出富集区域的区块，并以简单的 BED 格式报告合并后的坐标。或者，也可以报告 θc 以 bedGraph 格式保存每个区块，以保留读段密度估计的高分辨率细节。
4b. 点状读段分布：将后验均值后处理为峰候选区域
- 将背景速率定义为泊松分布 Pois(λa)，其值为所有读段计数（γ2）的平均值，并识别所有超过阈值 d 的区段。由于预期的点状峰具有更显著的富集，因此默认 δ 设为 Pois(λa) 的第 99百分位数。
- 将具有最大 θc 的区段设为候选峰顶点，并连接读段密度相似（允许±1个读段计数的微小变化）的侧翼区段。该连接区域被定义为候选结合位点。
- 将 ChIP 候选结合位点内的平均读段计数计算为 λ2，并针对输入背景进行假设检验，其中零假设 H0 为 λ1 ≥λ2，并根据 p 值阈值决定是否拒绝 H0 。以 BED 格式输出候选峰。
访问受限。请登录或开始试用以查看此内容。
结果
BCP 在识别组蛋白修饰数据中的广泛富集区域方面表现出色。作为参考，我们此前将我们的结果与现有工具 SICER3 的结果进行了比较，该工具已展现出较强的性能。为了充分展示 BCP 的优势，我们选择了一种已被深入研究的组蛋白修饰，以建立评估成功率的基础。基于此，我们分析了 H3K36me3，因其已被证实与活跃转录的基因体密切相关（图 1）。相比之下，H3K36me3 还被证明与抑制性标记 H3K27me3 呈互斥关系。我们进一步利用这些已知的关联关系，通过确定 BCP 所识别出的“岛区”与已知正相关和负相关区域的重叠比例，来展示其在识别准确性方面的优势。在此，我们通过更多高性能实例进一步证实 BCP 的优势。
我们先前的研究表明，BCP 检测出的富集岛尺寸（23.9 至 25.8 kb）明显大于 SICER 检测出的尺寸（2.7 至 10.7 kb）；较大的富集岛更符合 H3K36me3 富集区域通常呈广泛弥散型富集岛的传统预期 （PLoS Comp Bio, submi...
访问受限。请登录或开始试用以查看此内容。
讨论
我们着手开发一种用于分析ChIPseq数据的模型，该模型能够同样有效地识别点状和弥散状的数据结构。迄今为止，富集区域（特别是弥散区域）的识别一直较为困难，而这些区域反映了对较大岛屿尺寸的先验预期。为解决这些问题，我们采用了隐马尔可夫模型（HMM）技术的最新进展，该技术相较于现有的启发式模型以及创新性较弱的HMM具有诸多优势。
我们的模型采用具有显式公式的贝叶斯框架。这一特点与其他隐马尔可夫模型（HMM）有本质区别，使我们能够通过简单计算直接获得后验均值，即各片段的预期读段密度，而无需依赖耗时且计算成本高昂的模拟方法（如马尔可夫链蒙特卡洛方法）。因此，我们的计算时间和内存需求显著降低。在使用配备双核、2.0 GHz 节点及 2 GB 64 位内存的高性能计算集群分析约 2300 万条 H3K27me3 读段或约 2100 万条 H3K36me3 读段时，BCP 完成全基因组分析耗时不到一小时，而其他方法则需数小时至数天。这种时间节省仅需 2 GB 的 modest 内存即可实现。
此外，我们的模型对每个片段的各种均值进行条件...
访问受限。请登录或开始试用以查看此内容。
致谢
STARR 基金会奖项（MQZ），美国国立卫生研究院（NIH）资助项目 ES017166（MQZ），美国国家科学基金会（NSF）资助项目 DMS0906593（HX）。
访问受限。请登录或开始试用以查看此内容。
材料
| 姓名 | 公司 | 目录编号 | 评论 |
|---|---|---|---|
| 基于 Linux 的工作站 |
参考文献
- Park, P. J. ChIP-seq: advantages and challenges of a maturing technology. Nat. Rev. Genet. 10, 669-680 (2009).
- Barski, A., et al. High-resolution profiling of histone methylations in the human genome. Cell. 129, 823-837 (2007).
- Zhang, Y., et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 9, R137(2008).
- Zang, C., et al. A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics. 25, 1952-1958 (2009).
- Jothi, R., Cuddapah, S., Barski, A., Cui, K., Zhao, K. Genome-wide identification of in vivo protein-DNA binding sites from ChIP-Seq data. Nucleic Acids Res. 36, 5221-5231 (2008).
- Qin, Z. S., et al. HPeak: an HMM-based algorithm for defining read-enriched regions in ChIP-Seq data. BMC Bioinformatics. 11, 369(2010).
- Song, Q., Smith, A. D. Identifying dispersed epigenomic domains from ChIP-Seq data. Bioinformatics. 27, 870-871 (2011).
- Spyrou, C., Stark, R., Lynch, A. G., Tavaré, S. BayesPeak: Bayesian analysis of ChIP-seq data. BMC Bioinformatics. 10, 299(2009).
- Lai, T., Xing, H. A simple Bayesian approach to multiple change-points. Statistica Sinica. , (2011).
- Robertson, G., et al. Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing. Nat. Methods. 4, 651-657 (2007).
- Stitzel, M. L., et al. Global epigenomic analysis of primary human pancreatic islets provides insights into type 2 diabetes susceptibility loci. Cell Metab. 12, 443-455 (2010).
- Bernstein, B. E., et al. The NIH Roadmap Epigenomics Mapping Consortium. Nat. Biotechnol. 28, 1045-1048 (2010).
- Karolchik, D., et al. The UCSC Table Browser data retrieval tool. Nucleic Acids Res. 32, 493-496 (2004).
- Matys, V., et al. TRANSFAC: transcriptional regulation, from patterns to profiles. Nucleic Acids Res. 31, 374-378 (2003).
- Portales-Casamar, E., et al. JASPAR 2010: the greatly expanded open-access database of transcription factor binding profiles. Nucleic Acids Res. 38, D105-D110 (2010).
访问受限。请登录或开始试用以查看此内容。
重印与许可
标签
