ChIPpeakAnno：用于注释ChIP-seq和ChIP-ChIP数据的Bioconductor包

摘要

背景

ChIP和高通量测序（ChIP-seq）以及ChIP和基因组拼接阵列分析（ChIP-ChIP）已成为全基因组DNA结合蛋白靶点鉴定的标准高通量技术[1-4]。已经开发了许多算法和工具来分析ChIP-ChIP生成的大型数据集（参见[4])和ChIP-seq实验[1，5-10]。此类算法的输出通常以结合位点列表（也称为峰）的形式呈现，与对照样本相比，这些结合位点在ChIP样本中显著富集。所识别的结合位点通常被转换为一种格式，如BED或Wiggle（WIG），可以上传到UCSC基因组浏览器，这是一个开放访问、基于web的最新基因组序列数据源，与大量相关注释集成[11]。该资源允许用户构建自定义注释轨迹，以查看DNA结合位点与各种基因组特征（如基因、外显子、转录起始位点和保守元素）的接近程度。然而，搜索UCSC基因组浏览器对用户来说可能是一项艰巨的任务，尤其是当有大量结合位点或结合位点在基因组中广泛分布时。

已经开发了几个有用的web应用程序来管理和注释ChIP-ChIP数据[12-14]和ChIP序列数据[14]。然而，需要技术平台相关和基因组无关的批注释工具。这里我们描述了一个名为ChIPpeakAnno公司这有助于使用各种注释源对任何技术识别的结合位点进行批量注释，从而产生大量丰富的基因组区域，如ChIP芯片、ChIP-seq和CAGE。ChIPpeakAnno公司利用统计环境R/Bioconductor和各种注释源，如Ensembl、UCSC基因组数据库等。此外，用户可以灵活地用任何感兴趣的注释（如文献中的数据集）标记富集区域。此软件包可从Bioconductor获得，Bioconder是一个开源和开放开发软件项目，专门从事基于R（统计计算和图形系统）的生物数据分析和集成[15，16]。生物导体工具作为单独但可互操作的软件包分发，每个软件包专门用于生物数据分析的不同领域，例如利马微阵列数据分析软件包[17]以及生物反应器用于从联邦查询系统BioMart Ensembl检索基因组注释的包[11，18，19]。这个ChIPpeakAnno公司该软件包包含各种功能，用于批量通告从ChIP-seq、ChIP-ChIP或CAGE实验中识别的富集区域。

ChIPpeakAnno公司强调灵活性、集成性和易用性。用户可以从ChIP-seq、ChIP-ChIP、CAGE或任何实验中注释峰值，从而获得用户感兴趣的任何注释轨迹的染色体坐标列表。尽管其他软件程序中提供了一些功能，例如检索一组峰值的相邻序列，但ChIPpeakAnno公司在任何其他软件中都不可用。主要区别点CEAS，顺式基因组其他软件是ChIPpeakAnno公司允许将一组峰值与任何注释特征对象进行比较，例如与CpG孤岛进行比较，与保守元素（或未捕获的其他注释特征）进行比较欧洲原子能委员会http://ceas.cbi.pku.edu.cn/submit.htm或顺式基因组http://www.biostat.jhsph.edu/~hji/cisgenome/index.htm)（调查结果）或比较复制品之间的两组峰值（与NIH的Ivan Gregoretti进行个人交流）或复合物中的转录因子（未发表的数据）。此外，与ChIPpeakAnno公司，欧洲原子能委员会或顺式基因组未实施重叠显著性测试或基因本体（GO）富集测试。GO是一个描述基因产品的分子功能、生物过程和细胞分区的系统[20]。另一个主要优势是ChIPpeakAnno公司是插入其他注释包的能力/灵活性，例如生物反应器[17]和政府.db、ChIP-ChIP分析包，如林戈[21]和ACME公司[22]、其他快速发展的深层分析能力和基础设施（表​（表1）1)例如ShortRead（短阅读）[23],DEGseq公司[24]，边缘R[25],贝叶斯峰值[26],chipseq、ChIPseqR、Rolexa[27],BSgenome、IRanges、Biostrings、rtracklayer[28]，基因组图[29]和统计分析工具，如multtest和limma[17]生物导体（调查结果）。

可用性是ChIPpeakAnno公司。加载包后，一行代码(注释批中峰值)使用户能够从另一个数据集或任何感兴趣的注释轨迹中找到最近或重叠的特征，例如基因、外显子、miRNA、5'utr、3'utr和峰值。还为用户提供了动态获取注释的灵活性和功能(获取注释). 两行代码(获取丰富GO)允许用户查找丰富的基因本体术语。一行代码(制作维恩图)允许用户绘制维恩图并提供用于确定数据集之间重叠的重要性的p值。函数的重复调用查找重叠峰值使用户能够找到几个不同实验中的峰之间的重叠，这将有助于用户确定不同重复的峰如何重叠，以及不同转录因子的峰如何在复杂重叠中重叠。

实施

ChIPpeakAnno公司在R（统计计算和图形系统）中为ChIP-seq或ChIP-ChIP数据实现通用注释工作流[15，16]。为了促进组件重用和生物导体封装之间的兼容性，ChIPpeakAnno公司利用I范围包并将峰值列表表示为范围数据有效地找到最近或重叠的基因、外显子、5’utr、3’utr，microRNA（miRNA）或用户提供的其他自定义特征，如最保守的非编码元件、CpG岛或转录因子结合位点。所有调峰软件都会生成一个文件，其中至少包含一个染色体坐标列表ChIPpeakAnno公司包装需求。两张床http://genome.ucsc.edu/FAQ/FAQformat#format1和GFF（通用特征格式，http://genome.ucsc.edu/FAQformat#format3常见问题解答)是常见的文件格式，提供了将峰值或注释定义为数据行的灵活方法。因此，转换函数BED2范围数据和GFF2范围数据实现了将这些数据格式转换为范围数据对象。由于基因组注释定期/频繁更新，我们利用了生物反应器Bioconductor提供的软件包，可以从Ensembl中动态检索注释数据。为了快速访问，转录起始于常见基因组的位点，例如TSS.人类。NCBI36，TSS.人类。GRCh37，人类ExonPlusUtr。GRCh37、TSS.鼠标。NCBIM37，TSS.rat。RGSC3.4和TSS.斑马鱼。Zv8型作为预先构建的注释数据包包括在内。用户还可以灵活地从基因组特征包以及它们自己的注释数据，例如来自其他转录因子的结合位点列表、不同的ChIP制剂或不同的峰值调用算法。我们还利用了牛基因组包来实现允许检索与感兴趣的峰相关联的侧翼序列的功能。这有助于使用MEME等算法进行后续PCR扩增、克隆和/或模体发现[三，30]。为了确定确定的峰是否围绕具有特定GO项的基因富集，我们实施了GO富集测试。本测试采用超几何测试phyper公司在R中，并与来自政府.db包，特定于物种的GO注释包，例如组织健康状况数据库和多重性调整函数multtest公司生物导体封装。GO注释包每年更新两次。用批注的绑定站点ChIPpeakAnno公司可以导出为Excel文件，以便对较大的峰值列表进行轻松排序和统计分析。或者，相对于感兴趣的基因组特征（例如转录起始位点或外显子起始位点）的峰值分布可以很容易地绘制成饼图或直方图。此外，我们使用超几何检验实现了功能，以确定重复实验之间重叠的重要性、不同的峰值调用算法或复合体内转录因子之间的结合位点，并绘制维恩图以可视化重复之间重叠的程度。

结果

示例1：找到最近的基因以及到最近基因转录起始位点的距离

ChIP-seq或ChIP-ChIP分析的结果是一系列结合位点（作为染色体位置），与相应的对照样品相比，这些结合位点在ChIP样品中显著富集。下面的例子详细说明了如何在人类基因组中找到最近的基因以及与最近基因的转录起始位点（TSS）的距离，以获得一个结合位点列表（命名为我的峰值列表)类型为范围数据距离计算为结合位点的起始点与TSS之间的距离，即位于正向链上的基因的基因起始点和位于反向链上基因的基因结束点。

第一步是加载ChIPpeakAnno公司包、示例数据集和注释数据集。在本例中，示例数据集包含在未模拟单元格中识别的假定STAT1绑定区域[2]，注释数据集包含来自人类GRCh37的TSS坐标和链信息。

>库（ChIPpeakAnno）

>数据（myPeakList）

>数据（TSS.human.GRCh37）

在下一步中，函数annotatePeakInBatch（批内峰值）被调用以找到具有最近TSS的基因或不是最近TSS的重叠基因以及结合区域列表的相应距离。有时，峰值位于基因内，但远离基因的TSS。设置参数输出“both”输出TSS最近的基因和重叠的基因。参数最大间隙将最大间隙设置为重叠。参数倍数指示是否应为一个峰值返回多个重叠特征。

>annotatedPeak=annotatePeakInBatch（myPeakList，AnnotationData=TSS.human.GRCh37，output=“both”，multiple=F，maxgap=0）

带注释的峰可以保存为Excel文件，便于生物学家查看。

>write.table（作为.data.frame（annotatedPeak），file=“annotated PeakList.xls”，sep=“\t”，row.names=FALSE）

绘制峰值相对于TSS的分布图，可以鸟瞰峰值相对于感兴趣的基因组特征的分布。

>y=注释峰值$distancetoFeature[！is.na（注释峰值$distancetoFeature）&annotatedPeak$fromOverlappingOrNearest==“NearestStart”]

>hist（y，xlab=“到最近TSS的距离”，main=“”，breaks=1000，xlim=c（最小（y）-100，最大（y）+100））

从非刺激细胞中确定的假定STAT1结合区域生成的这种图([2])表明STAT1结合位点在TSS附近富集（图​（图1）。1). 还生成了如下饼图，以显示峰值相对于最近基因的相对位置分布（图​（图22).

在单独的窗口中打开

图1

人类STAT1结合位点相对于最近TSS的分布直方图是根据未刺激细胞中TSS周围20kb范围内的假定STAT1结合区域生成的[2]。该图显示STAT1结合位点在转录起始位点周围更加对称的区域富集。距最近TSS的平均距离为8533391±295725个碱基（平均值±SEM）。

在单独的窗口中打开

图2

人类STAT1结合位点相对于最近基因的饼图饼图由在未刺激的细胞中鉴定的推定STAT1结合区产生[2]。该图显示STAT1结合位点沿上游、下游和内部基因均匀分布。

>temp=as.data.frame（注释峰值）

>饼图（表（temp[作为.character（temp$fromOverlappingOrNearest）==“重叠”|（作为.charamer（temp$fromOverrappingOrNeorest）==%temp[以.character[作为.ccharacter（stemp$from OverlaappingOrnearest）==“重叠“，]$peak），]$insideFeature）中的“最近开始”&！temp$peak%

也可以使用获取注释调用之前的函数如下注释批中峰值。请参阅生物反应器可用列表的程序包文档生物艺术和数据集[18].

>mart=useMart（biomert=“ensembl”，dataset=“hsapiens_gene_ensembl“）

>注释=getAnnotation（mart，featureType=“TSS”）

为了用其他基因组特征注释峰值，有必要改变特征类型（例如，“exon”用于查找最近的外显子，“miRNA”用于查找最接近的miRNA，“5utr”用于查找最大的5'utr，“3utr”用来查找最大的3'utr）。也可以将自定义注释数据传递到函数中注释批量峰值例如，用户可能从文献、不同的生物复制、不同的峰值调用算法或作为转录复合物与研究中的蛋白质一起起作用的不同蛋白质中获得转录因子结合位点列表，并且有兴趣确定他/她的实验中峰值列表的重叠程度。调用函数之前注释批中峰值，有必要将两个数据集表示为范围数据，其中开始是绑定站点的开始，结束是绑定站点的末尾，姓名是绑定站点的名称，并且空间和搁浅分别是染色体名称和链，即结合位点所在的位置。

>myexp=RangedData（IRanges（开始=c（967654，2010897，2496704），结束=c（96 7754，2010997，24996804），名称=c（“Site1”，“Site2”，“站点3”），空格=c（”1“，”2“，”3“））

>文献=RangedData（IRanges（开始=c（9676592010898249670030758663123260），结束=c（96 78692011108249692030761663123470），名称=c（“t1”，“t2”，“t13”，“t4”，“tw”），空格=c（”1“，”2“，”3“，”1“、”2“），串=c（1，-1，-1,1）））

>annotatedPeak1=注释PeakInBatch（myexp，AnnotationData=文献，输出=“重叠”，倍数=F，maxgap=0）

峰值调用算法中文本格式的峰值可以很容易地导入到R中数据帧然后转换为范围数据。对于以BED或GFF格式表示的结合位点，BED2范围数据或GFF2范围数据用于将这些数据格式转换为范围数据.

示例2：确定重叠的重要性，并将重叠可视化为不同数据集之间的维恩图

第二个例子描述了如何确定重叠的重要性，在维恩图中可视化重叠，并从不同的数据集（例如不同的生物复制、不同的峰调用算法或作为转录复合物的不同蛋白质）中获得合并的峰。这里我们给出了使用不同生物复制品的示例。

第一步是加载ChIPpeakAnno公司包和三个示例数据集范围数据含有三个生物复制品酵母中确定的假定Ste12-结合区域[31].

>库（ChIPpeakAnno）

>数据（峰值.Ste12.副本1）

>数据（峰值.Ste12.Replicate2)

>数据（峰值.Ste12.Replicate3)

在下一步中，函数制作维恩图调用以生成维恩图，以可视化三个重复之间的重叠。此外，使用超几何检验进行了成对重叠显著性检验，并生成了p值。参数NameOfPeaks指示用于标记维恩图的数据集的名称。参数maxgap表示两个峰值范围之间的最大距离，以将其视为重叠。参数totalTest指示超几何测试中使用的潜在峰值总数.

>makeVennDiagram（RangedDataList（峰值.Ste12.复制1，峰值.Ste12.复制2，峰值.Ste12.复制3），NameOfPeaks=c（“复制1”，“复制2”，“复制3”），最大间隙=0，总测试=1580）

因此，生成了一个维恩图，用于可视化上述三个重复之间的重叠。成对重叠比较表明，从重复中确定的峰值重叠显著（图​（图3，三，p值<0.0001）。同样的分析也适用于Cse4的三个生物复制品，复制品1和复制品2之间的重叠在p值<0.01时显著，而其他两个重叠在p价值<0.05时显著（图​（图44).

在单独的窗口中打开

图3

酵母中三个生物复制体之间假定Ste12-结合位点的重叠维恩图是根据酵母中三个生物复制品的假定Ste12-结合位点生成的[31]。超几何检验表明，重复数据之间存在显著重叠（所有三个成对比较的p值<0.000001）。

在单独的窗口中打开

图4

酵母中三个生物复制体之间假定Cse4-结合位点的重叠维恩图是根据酵母中三个生物复制品的假定Cse4结合位点生成的[31]。超几何检验表明，生物复制品1和2在p值<0.001时重叠显著，而其他两个重叠在p值<0.05时显著。

复制品的峰值范围并不完全重叠。最好将重复数据中所有重叠的峰合并为合并的峰，覆盖重复数据中的所有重叠峰。调用函数查找重叠峰值可以生成合并的峰值。除了前面提到的参数外，还需要一个额外的参数倍数指示是否从多个重叠的峰值返回合并的峰值。

>合并峰值=查找重叠峰值（findOverlappingPeaks（Peaks.Ste12.Replicate1，Peaks.Ste12.Replicate 2，maxgap=0，multiple=F，NameOfPeaks1=“R1”，NameOf Peaks2=“R2”）$合并峰值，峰值。Ste12.Replicate3，maxgap=0，multiple=F，NameOfPeaks1=“Replicate1Repliate2”，NameOf Peaks2=“R3”）$合并峰值

接下来，通过使用SGD1.01注释合并的峰，可以获得最近的基因以及峰位置和最近的TSS之间的距离注释批中峰值示例1中所示的功能（图​（图55&6). 对Cse4结合位点进行了相同的分析（图​（图77&8). 非等方差t检验表明，Ste12-结合位点到最近TSS的距离分布（图​（图5，5，264±36个碱基）与Cse4结合位点（图​（图6，6，311±160个碱基）（p值=0.001）。Ste12-结合位点更多地分布在基因上游（图​（图55&6)而Cse4结合位点更多地分布在基因的内部和下游（图​（图77&8). 该结果与之前观察到的结果一致[31]。附加文件中提供了带注释的数据集1和附加文件2.

在单独的窗口中打开

图5

酵母中Ste12-结合位点相对于最近TSS的分布直方图是从酵母中确定的三个生物复制品合并的假定Ste12-结合区域生成的[31]。该图显示，Ste12-结合位点在转录起始位点上游和周围区域富集。到最近TSS的距离平均值为-264±36个基数（平均值±SEM）。

在单独的窗口中打开

图6

酵母中与最近基因相关的Ste12-结合位点饼图饼图是从酵母中鉴定的三个生物复制品合并的假定Ste12-结合区域生成的[31]该图显示，Ste12结合位点更多地分布在基因的上游和周围。

在单独的窗口中打开

图7

酵母中Cse4-结合位点相对于最近TSS的分布直方图是从酵母中发现的三个生物复制品合并的假定Cse4结合区域生成的[31]。该图显示Cse4结合位点在转录起始位点的内部和下游区域富集。到最近TSS的距离平均值为311±160个基点（平均值±SEM）。

在单独的窗口中打开

图8

酵母中Cse4-结合位点相对于最近基因的饼图饼图是从酵母中鉴定的三个生物复制品合并的假定Cse4结合区域生成的[31]。图中显示Cse4结合位点更多地分布在基因内部，并与基因末端重叠。

例3：获取结合位点周围的序列以进行PCR扩增或基序发现

第三个例子描述了如何获得结合位点周围的序列（在本例中，上游和下游序列为100 bp），用于PCR扩增、克隆或模体发现[三，30].

第一步是加载ChIPpeakAnno公司打包并创建示例peak数据集范围数据其次，特定于生物体牛基因组加载包，然后调用函数获取所有峰值序列.功能可用基因组显示特定于可用有机体的列表牛基因组数据包。在这个例子中，大肠杆菌数据包因其重量轻而被使用。

>库（ChIPpeakAnno）

>峰值=RangedData（I范围（开始=c（100500），结束=c（300600），名称=c（“peak1”，“peak2”）），空格=c（”{“类型”：“entrez-notide”，“属性”：{“文本”：“NC_008253”，“term_id”：“110640213”，“term_text”：“NC_008255”}}NC_008253号", "{“类型”：“entrez-notide”，“属性”：{“文本”：“NC_010468”，“term_id”：“170018061”，“term_text”：“NC_010468“}}NC_010468号"))

>文库（BSgenome.Ecoli.NCBI.20080805）

>peaksWithSequences=getAllPeakSequence（峰值，上游=100，下游=100，基因组=Ecoli）

要将序列转换为通用FASTA文件格式，调用以下函数。

>write2FASTA（peaksWithSequences，file=“test.fa”，width=50）

合并的Ste12结合位点的序列来自软件包英国基因组。Scerevisiae公司。UCSC.sacCer2（附加文件三). 通过运行MEME识别出重要的基序（E值<0.000001）[30]motif occurrence设置为ZOOP，最小宽度为8，最大宽度为20，其他参数为默认值（图​（图99).

在单独的窗口中打开

图9

酵母中Ste12-结合位点的基序.图案是使用MEME生成的[30]并且来自假定的Ste12-结合区的序列从酵母中鉴定的三个生物复制品中合并[31]。MEME的默认参数已选定，但motif occurrence设置为ZOOP，最小宽度设置为8，最大宽度设置为20。

结论

ChIPpeakAnno公司能够对从ChIP-seq、ChIP-ChIP、CAGE或任何技术中识别的结合位点进行批量注释，从而在统计编程环境R中使用现有注释数据为任何物种产生大量丰富的基因组区域。允许用户传递他们自己的注释数据，例如不同的ChIP准备和来自文献的数据集，或现有的注释包，例如基因组特征和牛基因组，提供了灵活性，同时与生物反应器包支持从BioMart数据库中检索最新注释。的主要优势ChIPpeakAnno公司能够/灵活地插入Bioconductor中的其他注释包、ChIP-ChIP分析包、其他快速移动的深层分析功能以及基础设施和统计分析工具。另一个优点是ChIPpeakAnno公司它使得能够在具有任何注释特征对象的一组峰之间、在来自复合物内的复制实验或转录因子的两组峰之间进行比较以及确定重叠的重要性。

这个ChIPpeakAnno公司该包以交互式手册的形式提供文档，说明各个函数的用法，以及包含演示面向案例的帮助会话的可执行代码片段的小插曲。vignette在包构建和安装时运行，因此也可以用作测试套件。本文中描述的一些示例也在小插曲中进行了演示。

可用性和要求

ChIPpeakAnno是GNU通用公共许可证v2.0下的一个开源软件包，并为生物导体项目做出了贡献。软件、源代码和文档可从下载http://www.bioconductor.org或通过键入source从R安装http://bioconductor.org/biocLite.R和生物灰岩（“ChIPpeakAnno”）。该软件包已经在OSX、Windows和各种Linux系统上进行了测试和运行。ChIPpeakAnno公司取决于R版本2.10。0或更高版本以及以下生物导体包：biomeRt、multtest、IRanges、limma、Biostring、BSgenome和GO.db此外，轻质有机体专用包装英国基因组。埃科利。NCBI.20080805号和组织健康状况数据库在构建期间安装，用于测试vignette中的代码片段。所有这些软件包都可以从Bioconductor下载，或使用http://bioconductor.org/biocLite.R脚本。

作者的贡献

LJZ起草了手稿。LJZ和HP开发了该软件包。CG、NDL、MRG、SML和DSL提供了科学建议。DSL执行MEME分析。所有作者都参与了写作并批准了最终稿。

补充材料

工具书类