一种用于多种ChIPseq数据类型全基因组分析的新型贝叶斯变点算法

以下实验的总体目标是利用染色质免疫沉淀测序数据中比对读段位置的密度，估算基因组范围内的后验平均读段密度。该目标通过预处理实现，即将比对后的ChIP-seq读段转换为分段密度谱，使相同数量的读段落入200个碱基对长度、互不重叠的分箱中。

作为第二步，将具有相同密度的相邻区间合并为更大的区块，然后在所有相邻区块的背景下，使用带有前向和后向滤波器的贝叶斯模型，递归计算每个区块的后验均值密度。其中，区块的读段计数由泊松分布建模，其 theta 参数服从伽马先验分布，该先验分布具有 alpha 和 beta 参数。接下来，通过判断每个区块的后验均值密度估计是否超过输入对照背景密度的第90百分位数，来评估其显著性，从而生成最终的富集基因组片段结果。这些结果展示了从原始测序读段到后验均值读段密度估计，再到 BCP 分析过程中 ChIP-seq 数据中富集“岛屿”的逐步演变过程。

此外，结果表明，BCP 的性能优于竞争对手工具 cer。与 CER 等现有方法相比，该技术的主要优势在于 BCP 采用了隐马尔可夫模型领域最新的 A 进展，因此相较于以往的启发式方法，能够更准确地刻画芯片数据解析中的细微差异。该方法有助于解决表观基因组学领域的关键问题，例如通过表征组蛋白修饰在全基因组范围内的富集模式，来揭示其生物学作用。

尽管这种严谨的方法可为 ChIP-seq 数据分析提供深入见解，但其基本框架也可应用于其他高通量测序数据的分析，例如在双硫测序（bisulfite sequencing）数据中鉴定差异甲基化区域、在 RNA-Seq 数据中识别新的转录位点、拷贝数变异分析，或多种微阵列平铺（microarray tiling）数据的分析。对该方法的可视化演示对于清晰理解其方法学原理及其优势至关重要。而理论上的优势则内嵌于软件之中。

此处演示的所有操作步骤均已整合到 BCP 软件包中的一个可执行程序中，该程序可在此视频中下载。文中描述了程序执行的步骤以指导软件运行。需要提供三个参数。

一个包含来自芯片样品的唯一比对读段的文件，以及一个相应的输入对照读段文件，还需提供一个输出文件名，以准备用于BCP分析的输入文件。首先，使用首选的短读段比对软件将测序运行产生的短读段比对至相应的参考基因组。比对到的位置应转换为六列的浏览器可扩展数据（BED）格式，每条比对读段对应一行以制表符分隔的数据，标明比对到的染色体、起始位置、终止位置、读段名称、得分以及链方向。

将芯片信号和输入图谱位置扩展至预设的片段长度。例如，DNA 经酶切或超声处理时所目标的片段大小，通常约为 200 个碱基对。随后，相邻区间内的片段计数将被合并汇总。

默认情况下，bin 大小设置为估计的片段长度，即 200 个碱基对。在一组具有相同读段计数的 bin 中，任何可能的变点最有可能出现在最外侧的边界处。因此，在两个读段计数相同的 bin 之间的内部边界上出现变点的可能性极低。

因此，将每条读段数相同的相邻区间归为一个区块。准备好输入文件后，只需输入屏幕底部显示的命令，即可调用BCP估计方法。每个区块的读段密度被建模为泊松分布，其均值参数theta服从伽马分布的混合分布，该混合分布具有alpha和beta参数，以及在任意区块处出现变点的先验概率。

以这种方式对每个区块的边界进行P条件处理，实际上相当于构建了一个无限状态隐马尔可夫模型（HMM）。超参数α、β和P通过最大后验似然法进行估计。对于每个区块θ_T的贝叶斯估计，显式计算其在给定Y_T条件下的期望值；传统上常用于HMM的前向-后向滤波方法虽然更为常见，但计算耗时，此处被更具计算效率的有界复杂度混合近似方法所取代，用于估计后验均值θ̂_T。所得后验均值将被平滑为近似的分段常数轮廓，因此具有相同θ̂_T值的区块应进一步合并，并更新其边界坐标。

BCP 使用每个区块的输入读段数量作为背景速率，并据此确定富集程度。通过一个简单的假设检验来判断染色质免疫沉淀位置在某个区块的平均密度是否超过某一显著性阈值。默认阈值为第90百分位数，在大多数情况下均适用。

BCP 随后会将相邻的、后验均值密度区块中超过富集阈值的部分合并为一个区域，并在浏览器中报告合并后的坐标。可扩展数据格式 BCP 在识别组蛋白修饰数据中的广泛富集区域方面表现出色。此处，BCP 的结果与 cser 的结果进行了比较；cser 是一种现有工具，在本实验室此前研究 H3K36 三甲基化的工作中表现出较强的性能，但相较于 BCP，其识别出的富集区域（“岛”）尺寸明显更小。

较大的岛状区域更符合传统的H3K36三甲基化富集区域广泛弥散的预期。然而，仅凭较大的岛状区域并不能单独说明准确性。因此，利用H3K36三甲基化岛状区域与活跃转录基因体的已知关联性，以及其与H3K27三甲基化岛状区域的互斥性，来评估BCP和CER的性能。与CER相比，BCP识别出更大的连续岛状区域，能更好地覆盖基因体，同时并未牺牲与H3K27三甲基化岛状区域重叠增加的程度。

BCP 通过 H3K36 三甲基化岛与活性基因保持高度重叠，其边界与基因体紧密对齐，且不会增加与转录受抑制的基因间区域基因或带有 H3K27 三甲基化抑制标记的基因之间的假阳性重叠程度。在评估 BCP 岛段识别在两组重复数据集中的可重复性时，观察到 BCP 并不像竞争算法 CER 那样严重依赖测序深度。BCP 具有稳健性和可重复性的进一步证据来自对其他不同区域的分析，即使测序深度降低，仍能展现出一致的岛段边界。为充分展示 BCP 的广泛适用性，获取了多种组蛋白修饰数据，包括点状分布的 H3K27 乙酰化、H3K9 乙酰化和 H3K4 三甲基化，以及弥散分布的 H3K9 三甲基化，此外还包括 H3K27 三甲基化和 H3K36 三甲基化。这些数据集均采用 BCP 和 CER 的默认参数设置进行分析。

在中心区域是PXDN基因上H3K36三甲基化的富集，标志着活跃的转录；在转录起始位点处，如预期般出现下降，同时还存在其他点状的活跃标志，包括H3K27乙酰化、H3K9乙酰化和H3K4三甲基化。PXDN基因下游紧邻的是由H3K27三甲基化富集标记的受抑制的基因间区域；而在相反侧的侧翼区域，则是一个被H3K27三甲基化标记的受抑制基因。再向外延伸一个区域。

我们的沉默染色质表现为存在 H3K9 三甲基化富集，这似乎表明 SN TG2 和 MYT1L 的沉默，其稳定性可能高于 H3K27 三甲基化所介导的抑制作用。该区域涵盖了在 ChIPseek 分析组蛋白修饰时所遇到的大部分现象。它展示了 BCP 的动态特性如何同时识别出点状的乙酰化修饰和 H3K4 三甲基化标记，并区分大片连续的 H3K27 三甲基化和 H3K9 三甲基化抑制区域，以及 H3K36 三甲基化所标志的活跃转录区域。

该算法的执行时间大约为30分钟，具体时长取决于读段数量和基因组特征结果。与其他方法通常所需的情况不同，本流程无需进行显著的优化步骤。利用BBCP可以研究染色质免疫沉淀的多种不同靶蛋白，包括各种组蛋白修饰以及DNA结合型转录因子，从而进一步探究表观基因组机制和基因调控相关的问题。

观看本视频后，您应该能够充分理解如何在 ChIP-seq 数据分析中利用 BCP 来识别组蛋白弥散修饰标记的富集区域。