内容

1简介

ChIP-seq作为检测蛋白质结合区域的标准技术已被广泛应用，其中专门开发了峰调用算法来服务于分析。由于测序技术，现有的峰值调用者缺乏对峰值重要性排序的权力，可能会产生序列上下文偏差，如。GC的偏见。gcapc旨在通过将GC效应建模为峰值调用来解决这一缺陷。gcapc还可以帮助细化其他峰值调用者所调用的峰值的显著性，或纠正一系列样本中预定义的基因组区域的读计数表的gc含量偏差。的gcapc包要求输入为一个ChIP-seq BAM文件(用于峰值调用/精炼)或一个读计数表(用于GC效果去除)以及其他可选参数。

对峰值调用/细化的常见分析包含四个步骤。

  1. 阅读报道。在此步骤中，BAM文件记录将分别转换为正向和反向链的基对分辨率的覆盖。

  2. 绑定宽度估计。该参数是ChIP实验中交联复合物中蛋白质结合区域大小的测量。同时，基于两股区域信号估计峰值检测的一半大小。

  3. GC效应估计。采用广义线性混合模型和EM算法来评估潜在的GC效应。

  4. 峰打电话/精炼。对于峰值调用，富集分数通过排列分析进行显著性评估。用浓缩分数和p值报告峰值。对于峰值细化，其他峰值调用者调用的峰值应该作为GRanges对象提供。新的富集显著值被添加为输入峰值的元列。

为了纠正计数表上的GC效应，基于函数的一步分析refineSites就足够了。

2开始

在R中加载包

库(gcapc)

3.准备输入

通过参考函数手册页，可以很容易地准备用于纠正计数表上GC效果的输入。在这里，我们主要关注用于峰值调用和细化的输入。输入可以最小到BAM文件的路径，这是用于排序读取的索引对齐记录。但是，鼓励指定其他选项以加速分析并提高准确性。以下是用户可以自定义的选项集。

  1. BAM记录过滤选项。在函数中read5endCoverage，读取可以过滤所选染色体，映射质量，重复去除等。如果只分析染色体的一个子集，则可以大大加速下游分析。如果一个ChIP-seq实验被深度测序，这实际上表明了一种分而治之的策略。在这种情况下，基于每个染色体级别的分析可以节省大量内存。

  2. 片段排序选项。如果对测序片段的大小有先验知识。函数中的可选参数bindWidth可以指定将搜索范围限制在较窄的范围内;或者，如果预先知道绑定宽度，则可以省略此函数。请注意，这个结合宽度可能并不等同于生物学中蛋白质的结合宽度，因为它可能受到交联操作的影响。

  3. 用于GC效应估计的样本量。默认值为0.05，这意味着如果基于全基因组进行分析，将使用5%的基因组。然而，对于较小的基因组或较小的染色体子集，这个大小应该调高以确保准确性。或者，对基因组进行多次采样，并使用平均估计来避免抽样偏差。注:样本量越大，采样次数越多，GC效应估计的计算时间越长。

  4. GC影响估计的GC范围。如手册页所示，GC范围为(gcrange参数)应仔细选择。原因是GC含量极低或极高的区域有时会作为异常值，当混合模型拟合中选择的前景区域过少时，可以驱动回归线。当所研究的结合蛋白的全基因组结合事件太少时，就会发生这种情况。

  5. EM算法的先验和收敛性。可以调整EM算法的选项以加速迭代。

  6. 排列。正如我们在函数帮助页面中所建议的，适当的排列时间可以节省时间，并确保准确性。

在这个小插图中，我们将使用嵌入式文件chipseq.bam作为一个例子来说明这个包。该文件包含大约80000个来自人类21号染色体CTCF ChIP-seq数据的读取。

Bam <- system。文件(“extdata”、“chipseq。bam”,包= " gcapc”)

4峰打电话/精炼

具体算法请参考我们的论文(Teng and Irizarry 2017)．

4．1读取范围

第一步是为正向和反向链生成读取覆盖。覆盖范围基于单核苷酸分辨率，仅使用BAM记录的5 '端。这意味着，如果不允许重复，每个核苷酸的最大覆盖范围是1。

cov <- read5endCoverage(bam) cov ## $ fwd# # RleList of length 1 ## $chr21 ## integer-Rle of length 48129895 with 40225 runs ##长度:9414767 1 8350 1…1 116 1 41437 ##值:0 1 0 1…1 0 1 0 ## ## ## $rev ## RleList长度为1 ## $chr21 ##整数-长度为48129895的rle与40427运行##长度:9412972 1 3087 1…1 367 1 34767 ##值:0 1 0 1…1 0 1 0

Obejct浸是两个元素的列表，分别表示正向链和反向链的覆盖范围，而每个元素是单个染色体上的覆盖范围。

4.2绑定的宽度

第二步是估计ChIP-seq实验的绑定宽度和峰值检测半窗大小。如果事先知道绑定宽度，则可以省略此步骤。进一步将绑定宽度作为区域单元的大小，用于有效的GC偏差估计和峰值调用。峰值检测半窗大小用于定义侧翼区域的宽度。

如果从测序片段中获得更多信息，则可以加快这一步骤。例如，缩小范围大小会有所帮助。

bdw <- bindWidth(cov, range=c(50L,300L)， step=10L) ##开始估计bdwidth。# #……循环1绑定宽度估计## ......循环2绑定宽度估计## ......循环3绑定宽度估计## ......循环最终绑定宽度估计## ......估计绑定宽度为111 ## ......估计峰值窗口一半大小为220 ##从两个股细化峰值窗口一半大小的区域## ......................# #……细化峰值窗口一半大小为135 bdw# # [1] 111 135

4.３GC效应

该步骤使用所提出的模型执行GC效应估计。值得注意的是，通过允许显示图表，人们可以查看中间结果，这些结果为您提供了对ChIP-seq数据的直接感觉，例如GC影响的程度。此外，默认情况下启用EM算法迭代，以显示日志可能性变化的跟踪，并打印其他通知消息。

layout(矩阵(1:2,1,2))gcb <- gcEffects(cov, bdw, sampling=c(0.25,1)， plot=TRUE, model='poisson') ##开始估计GC效果## ......'supervise'中的范围太少/太多。选择随机抽样## .........对25%的地区进行1次抽样，共108395个地区## ......计数是## ......计算GC内容与侧翼79 ## ......估计GC效果## .........迭代1 ll -35148.61增量109906.9 ## .........迭代2 ll -34316.64增量831.9672 ## .........迭代3 ll -34285.41增量31.23412

在这里，25%的窗口采样1次重复。左图使用估计的绑定宽度作为区域单位，提供了正向和反向链信号之间的相关性。右图显示了使用混合模型的原始和预测GC效应。

这里需要注意的其他两个选项是监督而且模型．如果监督选项被指定为GRanges对象，它提供了一组潜在的峰值，并允许更有效的采样过程。具体来说，这两种混合物分别从前景(信号)和背景区域采样。模型选项允许在模型拟合中在泊松分布和负二项分布(默认)之间切换。理论上，负二项式假设比泊松假设更准确。尽管如此，泊松是一个很好的近似于负二项式的GC效应估计，并显示出比负二项式更快的计算速度，特别是当选择的箱子总数很大时。

4.4最大的意义

这是高峰呼叫的最后一步。它使用在前面步骤中生成的信息，计算富集分数并执行排列分析来提出重要的峰值区域。最终形成的山峰农庄对象，元列用于记录重要性。还会打印其他通知消息。

layout(矩阵(1:2,1,2))peaks <- gcapcPeaks(cov, gcb, bdw, plot=TRUE, permute=100L) peaks <- gcapcPeaks(cov, gcb, bdw, plot=TRUE, permute=50L)

峰值与413范围和2 # #农庄组织对象元数据列:# # seqnames范围链| es pv # # < Rle > < IRanges > < Rle > | <数字> <数字> # # [1]chr21 9827175 - 9827175 * 23.063 | 9.17473 e-08 # # [2] chr21 9909602 - 9909602 * 8.042 | 1.65546 e-02 # # [3] chr21 9915442 - 9915442 * 6.937 | 2.80852 e-02 # # [4] chr21 11175039 - 11175039 * 6.978 | 2.75415 e-02 # # [5] chr21 15626034 - 15626034 * 32.631 | 9.17473 e-08  ## ... ... ... ... . ... ...[409] chr21 47705703-47705851 * | 6.330 3.72794e-02 ## [410] chr21 47745293-47745459 * | 12.006 1.83045e-03 ## [411] chr21 48055040-48055206 * | 13.856 5.48833e-04 ## [412] chr21 48059797-48059983 * | 11.168 3.03785e-03 ## [413] chr21 48081170-48081315 * | 33.969 9.17473e-08 ## ------- # seqinfo: 1个来自未指定基因组的序列;没有seqlengths

值得注意的是，这里使用不同排列次数的两次测试的浓缩分数的截止值几乎相同，这表明允许少量排列以节省时间。左图显示的是基于50次排列的富集分数的截止值，右图显示的是基于100次排列的富集分数的截止值。请注意，我们仅使用21号染色体进行说明，因此将置换次数从默认的5增加到这里的50。

4．5峰精炼

为了消除GC对其他峰值调用者输出的影响，这个包提供了针对给定峰值细化富集意义的函数。峰值必须作为GRanges对象提供。为了减少潜在的假阴性，建议采用一组灵活的峰，这意味着建议包括显著峰(例如p<=0.05)和不显著峰(例如p>0.05)。如果峰值总数不太大，则合理的峰值集包括其他峰值呼叫者p-value/FDR小于0.99的所有峰值。

newpeaks <- refinePeaks(cov, gcb, bdw, peaks=peaks, permute=50L) plot(newpeaks$es,newpeaks$newes,xlab='旧分数'，ylab='新分数')

newpeaks # #农庄与413范围和4元数据对象列:# # seqnames范围链| es光伏新# # < Rle > < IRanges > < Rle > | <数字> <数字> <数字> # # [1]chr21 9827175 - 9827175 * 23.063 | 23.063 - 9.17473 e-08 # # [2] chr21 9909602 - 9909602 * 8.042 | 8.042 - 1.65546 e-02 # # [3] chr21 9915442 - 9915442 * 6.937 | 6.937 - 2.80852 e-02 # # [4] chr21 11175039 - 11175039 * 6.978 | 6.978 - 2.75415 e-02 # # [5] chr21 15626034 - 15626034 * 32.631 | 9.17473 e-08 32.631  ## ... ... ... ... . ... ... ...## [409] chr21 47705703- 47745459 * | 12.006 1.83045e-03 12.006 ## [411] chr21 48055040-48055206 * | 13.856 5.48833e-04 13.856 ## [412] chr21 48059797-48059983 * | 11.168 3.03785e-03 11.168 ## [413] chr21 480597170 -48081315 * | 33.969 9.17473e-08 33.969 ##新pv ## <数字> ## > ## >.19522e -02 ## [5] 1.23901e-07 ## # ... ...## [409] 6.53721e-02 ## [410] 5.49960e-03 ## [411] 2.09207e-03 ## [412] 8.32764e-03 ## [413] 1.23901e-07 ## ------- # seqinfo:来自未指定基因组的1个序列;没有seqlengths

在这里，两个新的元列被添加到之前的峰值(GRanges)中，包括调整后的显著性和p值。值得注意的是，新的富集分数实际上与以前计算的相同(上图)，因为峰值区域以前是由gcapc．在实际应用中，其他调峰者所调用的峰值区域及其意义大多不同gcapc是否有较强的GC偏倚。如果细化这些峰值显著性，改善应该是明显的，正如我们在论文中所展示的那样。

5校正计数表的GC影响

在这个包中，函数refineSites如果有些人更感兴趣的是校正预定义区域的GC效应，而不是峰值调用。在我们的论文中，我们使用这个函数来调整ENCODE报告站点的信号(Teng and Irizarry 2017)．该函数易于使用，计数表和相应的基因组区域是两个必需的输入。有关使用此函数的详细信息，请阅读函数手册页。

6总结

在这篇短文中，我们介绍了这个包中的主要功能，并说明了它们是如何工作的。通过遵循这些步骤，用户将能够消除对峰值识别或读计数信号的潜在GC影响。

参考文献

滕明翔，Rafael A. Irizarry, 2017。“计算gc含量偏差可减少芯片序列数据中的系统误差和批量效应。”基因组研究．